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ABSTRACT 

The detection and characterization of self-organized criticality (SOC), in both real and simulated data, 
has undergone many significant revisions over the past 25 years. The explosive advances in the many nu¬ 
merical methods available for detecting, discriminating, and ultimately testing, SOC have played a critical 
role in developing our understanding of how systems experience and exhibit SOC. In this article, methods 
of detecting SOC are reviewed; from correlations to complexity to critical quantities. A description of 
the basic autocorrelation method leads into a detailed analysis of application-oriented methods developed 
in the last 25 years. In the second half of this manuscript space-based, time-based and spatial-temporal 
methods are reviewed and the prevalence of power laws in nature is described, with an emphasis on event 
detection and characterization. The search for numerical methods to clearly and unambiguously detect 
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SOC in data often leads us outside the comfort zone of our own disciplines - the answers to these ques¬ 
tions are often obtained by studying the advances made in other fields of study. In addition, numerical 
detection methods often provide the optimum link between simulations and experiments in scientific re¬ 
search. We seek to explore this boundary where the rubber meets the road, to review this expanding field 
of research of numerical detection of SOC systems over the past 25 years, and to iterate forwards so as to 
provide some foresight and guidance into developing breakthroughs in this subject over the next quarter 
of a century. 

Subject headings: Self Organized Criticality, numerical methods 


1. INTRODUCTION 

Self-Organized Criticality (SOC) is a statistical property of many time-varying systems. Aschwanden et al. 
(2014) (this volume of SSR) present a detailed description of SOC in solar and astrophysical settings; for the purposes 
of this current paper, SOC is considered in the wider aspect of any physical system that displays the scale invariance 
in both time and space leading to a critical point. It is often observed in slowly driven, but non-equilibrium, systems 
and, perhaps most importantly, complexity naturally arises in the system without any fine-tuned parameters as input. 
Although well-known earlier work {e.g., Neumann 1966; Mandelbrot 1975) had shown that complexity could arise 
from simply-governed, slowly driven systems, the seminal paper of Bak et al. (1987) provided the breakthrough in 
this subject by showing that all the so-called SOC features {e.g., fractal geometry, scale-invariance, power laws) arise 
from simple systems and lead to a critical point with no fine tuning of the input. Hence the system is both self- 
organized and critical. The large volume of research resulting from Bak et al. (1987) includes many articles on how to 
recognize SOC in a system. It is the 25 years of these numerical detection methods that we review in this paper. 

The power of SOC lies in the ability to both describe and explain a large variety of physical systems in a quan¬ 
titative and physically-motivated manner . From sand piles (Bak et al. 1987) to solar flares (Lu and Hamilton 1991), 
from fractures (Turcotte et al. 1985) to forest fires (Drossel and Schwabl 1992); from asteroids (Ivezic et al. 2001) 
to accretion disks (Dendy et al. 1998), SOC provides a mathematically tractable and understandable route to study 
complex systems. The scale-free, dimensionless, nature of SOC conveniently encompasses much of the universe. 
The concept of simple beginnings - assuming a starting grid and apply a few rules regarding distribution of excess 
amongst nearest neighbors - is an attractive model to many scientists, spanning subjects from physics and chemistry 
to economics and sociology. However, every SOC researcher ultimately reverts back to the same set of unanswered 
questions - How can I tell whether my system is truly SOC, or if it is just displaying SOC-like behavior? How can 
I detect SOC in such a way that I can confidently distinguish it from other potential physical sources? The route to 
answering these questions begins in Section 2.1 with the seemingly-simple studies of autocorrelations, described in 
terms of symmetries leading to diffusion models, and correlations functions leading to surface growth models. We 
end this discussion with a detailed look at the methods of measuring correlation functions, with a emphasis on the 
Manna model. The models introduced in this section are all guided by simple sets of rules of particle interaction gov¬ 
erning how particles spread apart {i.e., diffusion), how particles clump together {i.e., growth), and the redistribution of 
particles upon reaching a threshold value. In Section 2.2 we move from a discussion of products of field values {i.e., 
correlation functions) to a discussion of increments {i.e., structure functions). The value of the structure function as a 
complementary approach is highlighted with respect to determining linear ranges in log-log plots, with an application 
to solar magnetic fields. Application-oriented methods (Section 2.3) provide a third approach to numerical detection 
of SOC. We end our discussion of numerical methods in Section 2, by studying the advantages of block-scaling as a 
sub-sampling method to be used when little data is available to the scientist. 
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With this toolkit in hand. Section 3 contains a review of the many approaches developed over the last 25 years 
to identify individual SOC features and events. We split these studies into the three areas of spatial, temporal, and 
combined spatio-temporal. By performing this three-way split we merely seek a convenient route to provide some 
narrative to the reader; we do not suggest that these techniques differ in some fundamental way. When studying 
images in Section 3.1 we usually require thresholding, and considerations of 3D volume. As we typically only have 
2D images, this consideration leads us to discuss the potential 2D fingerprint of a 3D SOC system. As a follow-on 
from this type of thought process, one need only look at that most common feature of SOC detections of power laws 
in Section 3.2. It is clearly trivial to plot data on a log-log set of axis and find a straight line fit. The real purpose 
of this scientific endeavor should be research performing a set of logical deductive steps showing that such data are 
truly described by a power law, and that this power law can only be the result of an SOC system. The discussion in 
Section 3.2.1 shows how rarely we achieve such a scientific nirvana. Only when we fully comprehend issues such as 
the detection of power laws, and issues of data sampling and pulse pile-up can we then move to discuss waiting-time 
distributions as a possible signature of SOC. We conclude in Section 3.3 by showing how spreading and avalanche 
exponents provide vital tools to study spatio-temporal structures, with an emphasis on examples from magnetospheric 
and solar physics. 
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2. Methods of numerical detections of SOC 

The basic approach to test for the existence of SOC in numerical or observational data is to extract a series of 
events and test if these features are in some way connected. Events can are often called features, clusters, storms, 
objects, explosions, instabilities - the nomenclature is often different but the principle is the same. In Section 3 we 
will proceed to perform a synthesis on methods of extraction of these events, however here in Section 2 we first review 
existing methods of testing for connections between events, starting with the autocorrelation function and its modern 
extensions (Section 2.1), moving onto structure functions (Section 2.2) and then focusing on application-oriented 
methods developed in the last 25 years (Section 2.3). 


2.1. Autocorrelation functions 

Autocorrelation functions have a long history in the study of critical systems (Stanley 1971). While they are 
defined on the microscopic scale, they bridge the gap to the large scale and typically display scaling on these larger 
scales in both space and time. As such, correlation functions are at the heart of the theoretical description of scaling 
phenomena in systems with many interacting degrees of freedom, yet numerically and experimentally they are often 
inaccessible. The provision of numerical detection methods for the study of SOC systems hinges critically on a 
fundamental understanding of correlations functions in the study of traditional systems. In the following section, 
correlation functions are introduced in broad terms, highlighting some basic features and symmetries that are important 
for a later discussion of SOC. Readers familiar with these two topics may wish to skip to Section 2.1.3 where we 
discuss some basic null models in order to motivate the focus on some characteristics of correlations often found in 
non-trivial systems exhibiting SOC. Some parallels are drawn from the study of surface growth and interfaces and then 
the basic measurement methods are exemplified using the Manna Model (Manna 1991). 

SOC systems evolve in time and extend in space due to the interaction of their local degrees of freedom (local 
activity of avalanching, energy, particle density, height etc.). The propagation of this interaction in time and space 
can be captured by autocorrelation functions. As SOC systems demand evolution to a critical point, it is expected that 
every part of a system interacts with every other part of a system, as well as with their history, in such a way that does 
not allow for degrees of freedom to be dropped on the basis that they are too remote in space or time. Even the most 
local features cannot be studied in an isolated fashion, as local degrees of freedom self-interact, mediated by their 
environment. Correlation functions are therefore used to both measure and quantify these effective interactions at the 
most basic level. 


2.1.1. Basic features 

The most basic autocorrelation function of local degrees of freedom ^(r,t), such as the local particle density, 
energy, magnetization etc ., at position r and time t is 

C(r2,f2,ri,ti) = (^(r2,f2)^(ri,fi))-(0(r2,f2))(^(ri,fi)) (1) 

where (•) takes the expectation value, i.e., it is the ensemble average. If (j){r 2 f 2 ) 0(riTi) uncorrelated, in 
particular when they are independent, the joint probability density of (j){r 2 ,t 2 ) and 0(ri,ti) factorizes and therefore 
(0(r2,f2)0(ri,fi)) = (0(r2,f2)) (^(ri,fi)) (i.e., the correlation function vanishes) C(r2,f2,ri,fi) = 0. This is obvi¬ 
ously a rather trivial situation - correlations do not matter for these types of degrees of freedom and the behavior of 
one is not influenced by the behavior of any other. When ri = r 2 and fi = t 2 the correlation function C(r 2 , f 2 , ri Ti) in 
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fact describes the variance of the local (j). Alternatively C{r 2 , 127 ^ 1 ,h) may be thought of as a measure of fluctuations 
relative to the background as Eq. (1) can be re-written as 

C(r2,f2,ri,fi) = ^(0(r2,f2)-(0(r2,f2)))(0(ri,fi)-(0(ri,ti)))^ . (2) 

The result is large when large fluctuations at ri,fi match large fluctuations at r 2 ,f 2 , and it is small when they typically 
miss each other. The correlation function might be negative, signalling anti-correlations if positive fluctuations at ri, t\ 
typically occur when they are negative, ^(r 2 ,f 2 ) ~ ( 0 (r 27 f 2 )) < 0, at r 2 ,t 2 . 


2.1.2. Symmetries 

Symmetries may simplify the dependence of C{r 2 ,t 2 ,ri,ti) on the two points in both space and time. If the 
system is translationally invariant, thenC(r 2 ,f 2 ,ri,fi) is afunction only of the difference r 2 —ri, i.e., C{r 2 ,t 2 ,ri,ti) = 
C(r 2 — ri,f 2 , 0 ,ti). If it is, in addition, invariant under rotations, then it is only a function of the distance |r 2 — ri|. 
When estimating C(r 2 ,f 2 ,ri,ti) from numerical or observational data, these invariances can be used to improve the 
estimates, for example in the form 

C'(r,f 2 ,fi) /d^r'C(r',r' + r,f 2 ,fi) (3) 

Jv 

where the integration runs over the entire t/-dimensional volume V of the system. A system with boundaries cannot 
be expected to be truly translational invariant, so this is often used as a suitable approximation only in relatively 
small localizations deep inside the system. Most SOC systems require boundaries in order to dissipate energy or 
particles driven into it, and they are often not translational or rotational invariant, although some basic symmetries, 
{e.g., due to the shape of the system) remain. A typical example is an inversion symmetry about the origin, so that 
C(r 2 ,t 2 ,ri,fi) =C(-r 2 ,f 2 ,-ri,ti). 

Similar simplifications apply in the time domain. If correlation functions are translationally invariant in time 
the system is said to be stationary, i.e., C(r 2 ,f 2 ,ri,fi) = C(r 2 ,t 2 — ti,ri,0) = C(r 2 , 0 ,ri,fi — f 2 )- By construction of 
Eq. (1), C is invariant under permutations of the indices, C(r 2 ,f 27 ri ;fi) = C(ri ,r 2 ,ti,t 2 ). If C is additionally invariant 
under rotation and translation, C(r 2 ,f 27 ri 7 li) = C’(ri,f 2 ,r 2 ,ti), then, by definition, Eq. (1) implies invariance under a 
change of sign of f 2 — ti, 

C(r2,f2-fi7ri,0) =C(r 2 ,t 2 ,ri,ti) =C(ri,ti,r 2 ,f 2 ) =C(r 2 ,ti,ri,f 2 ) =C(r 2 ,fi -t2,ri,0) . (4) 

However, correlation functions are often of the form 

G(r 2 ,f 2 ,ri,ti) = (^(r 2 ,f 2 )r(ri,fi))-(^(r 2 ,t 2 ))(r(ri,fi)) (5) 

where v/(ri,fi) denotes a perturbation of the system at time t\ and position ri and (j){r 2 ,t 2 ) is the response at time t 2 
and position r 2 . In this case a change in the sign of f 2 — f 1 reverses the causal order and therefore the correlation function 
G(r 2 ,f 2 ,ri,fi) is not invariant under that change, as it is not invariant under an exchange of indices - G(r 2 ,t 2 ,ri,ti) 
G(ri,fi,r 2 ,f 2 ), as they refer to different entities. Initial conditions generally play the same role as perturbations or 
boundary conditions - the presence of initial conditions undermines stationarity and time reversal symmetry, just as 
the presence of boundary conditions undermines translational invariance and inversion across arbitrary points. In 
order to distinguish Eq. (1) from Eq. (5) in the context of SOC, the former is often referred to as the activity-activity 
autocorrelation function and the latter, less common, is referred to as the propagator or response (correlation) function. 
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Two-point correlation functions are simply correlation functions evaluated at two sets of coordinates (or, if suit¬ 
able symmetries are found, differences of two sets of coordinates). In most applications, two-point correlation func¬ 
tions are either evaluated at the same time fi = t 2 , known as equal time correlation functions, or at the same point in 
space ri = r 2 , and known as temporal correlation functions or two-time correlation function. The behavior captured 
by an equal time correlation function is thought to be due to a common source, like the simultaneous ripples on the 
surface of a pond at two points are caused by a stone dropped at the origin {e.g., it is very instructive to study correla¬ 
tions in a deterministic system as simple as 0(r,t) = sin(ko|r| — (Bot)/|r| for some fixed ko and (Oo). If the correlation 
function is intended to measure causal relationships, such as in Eq. (5), it must necessarily vanish at equal times for 
rj ^ r 2 , as a perturbation is expected to require time to propagate from ri to r 2 . To stay in the same picture, the 
response function in Eq. (5) would measure the response at r 2 ,f 2 to a stone dropped at 


2.1.3. Basic diffusion examples and null models 


In many cases, the field ^ denotes a particle density and the null-models of correlations in time and space are 
Poisson and Gaussian processes. The former refers to processes where events occur completely independently with 
constant rate, the latter to the random and interaction-free spreading of a quantity subject to conservation and conti¬ 
nuity. In the former case, all connected correlation functions vanish. In the latter case, plain diffusion with constant 
Brownian diffusion coefficient D introduces correlations between different points in time and space. If a single, freely- 
diffusing particle is created at time fi and position ri, the relevant correlation function in d Euclidian dimensions is 
(van Kampen 1992; Strauss 2007), 

/ \ j.j2 

G(r2,t2,ri,fi) = 0(f2-Ii) ( , -, ) e . (6) 


It describes the expected particle density at r 2 ,f 2 following the creation of the particle at ri,fi. Equivalently, it is the 
probability density of finding that particle at r 2 ,f 2 after it has been created at ri,ti. Eq. ( 6 ) is also the solution of the 
deterministic diffusion equation. 

Hwa and Kardar (1989) proposed a model more relevant to SOC by introducing a source 77 (r,t), so that 0 (r,f) = 
Jd'^r' /Qdf'G(r,f,r',t')77(r',fO- V describes Gaussian white noise with some amplitude 2r^, then is the 

height of an interface subject to Edwards-Wilkinson dynamics (Edwards and Wilkinson 1982; Krug 1997). It can be 
thought of as a surface, or a diffusive field, relaxing under the influence of surface tension v =D, while being exposed 
to random addition and removal of material (parameterized by E^). In one dimension the equal time correlation 
function becomes 


C(r2,ri,t,t) = 2r2 





|r2-ri 



|r2-ri| 


(7) 


and the temporal correlation function starting from a flat interface is then 


C(r,r,f,t) = 2r^ 


yti+ti — a/ig— fii 


( 8 ) 


In terms of observables, this is what is typically studied in SOC systems - namely the correlation of the local height or 
the particle numbers between sites. 

It is important to note that the distinction between the response function, G, and the correlation function, C, is 
more than a technicality. The former is the correlation function for the propagation of a perturbation within the degrees 
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of freedom - it addresses the question of how the degrees of freedom, the field <j), reacts to a perturbation. The latter, 
on the other hand, describes the correlations seen in the degrees of freedom as the system evolves. These are mediated 
by the propagator that communicates events, in particular any external driving, to other sites in the system. To draw 
a rough parallel to seismic events: G is the seismic signal measured r 2 ,f 2 throughout the Earth’s crust as a bomb 
detonates at ri, fi, whereas C are the correlations between the signal at r 2 , f 2 and ri, fi as the earth crust evolves under 
its natural dynamics. 


2.1.4. Temporal and spatial correlations 

Long-range temporal correlations are frequently found in non-equilibrium systems, even when the microscopic 
interaction is trivial in the technical sense discussed below (Grinstein 1995). Even directed models display scaling 
in temporal correlation functions (Pruessner 2004b). Non-trivial spatial, correlations are generally regarded as the 
signature of interactions that dominates the large scale. Temporal correlations are often quantified by the correlation 
time T (see also the correlation length ^ introduced below). The correlation time is defined by the asymptotic decay 
of the correlation function C(r,t 2 ,r,fi) exp(—|f 2 — L for large \t 2 — t\\. It can be defined in a correspondingly 
similar fashion for the propagator, or response function, G(r,f 2 ,r,fi) exp(—|t 2 — L/t)). This structure follows 
necessarily if the observable ^(r,f) is subject to Markovian dynamics, so that T is in fact determined by the negative 
inverse logarithm of the second largest eigenvalue of the Markov matrix (van Kampen 1992). 

An equation very similar to the Edwards-Wilkinson equation was suggested by Hwa and Kardar (1989) as a 
description of SOC phenomena with a possible mass term, e, that parameterizes an attenuation of the signal. The 
resulting equal-time correlation functions \nd = l and c/ = 3 dimensions are, in the limit of large times, 

limCi(r2,f,ri,0 = (9a) 

limC3(r2,f,ri,f) = ^ (9b) 

2v|r2 —ri| 

These are also known as Ornstein-Zernike-type correlation functions - namely Eourier transforms of T^/(vk^ H- e) 
obtained in the Ornstein-Zernike approximation (Stanley 1971, Chap. 7.4.2, Barrat and Hansen 2003, Chap. 5) for 
the structure function in liquids. In some settings, studying the Eourier transform in space, essentially produces the 
structure factor whereas studying the Eourier transform in time, essentially produces the power spectrum (Abramenko 
et al. 2003). 

The examples above are instances of trivial correlations in different disguises. Apart from the fact that the 
only scale mentioned is that of the diffusion constant D or the surface tension v, which imposes the typical relation 
between time and space f r^, it is the triviality in the technical sense that makes them proper null-models. Trivial here 
means that the correlations are produced in the absence of interaction, which, in turn, is absent because the processes 
considered above are linear, i.e., the stochastic partial differential equations of motion are linear in the field (j). The 
equivalence of linearity and lack of interaction can be understood by noticing that solutions can be superimposed - 
adding one solution to another produces a new solution. In other words, the solution to an initial condition with two 
particles initially deposited is just the sum of the solutions for each particle individually - the particles do not see 
each other. Therein lies the reason for the interest of statistical mechanics in non-trivial, spatial correlations. Their 
space-dependence is normally quantified by matching correlation functions to the scaling form. 


( 10 ) 



- 8 - 

with a so-called metric factor a (independent of ^, see Christensen et al. 2008), Euclidean dimension d, universal ex¬ 
ponent 77 , also known as the anomalous dimension, and a scaling, or cutoff, function that describes how correlations 
eventually decay on a scale beyond the correlation length, ^. The divergence of the correlation length at the critical 
point is probably the most direct signal of criticality. In SOC, where systems are expected to organize themselves to 
the critical point, the correlation length is naturally limited by the system size L and all scaling of global, system-wide 
observables in SOC is therefore finite size scaling (Barber 1983). As such, one of the most direct tests of the system 
being at criticality is to demonstrate that ^ L. 

Eq. (10) is not normally expected to hold on short scales, where lattice effects become important. Rather, it 
describes an asymptotic behavior in large distances |r| and for large correlation lengths In particular, it is not 
expected to capture the degeneration of C(r,f,0,t) into the variance at r = 0. Even when the exponent becomes 
negative, —(d — 2 + ri) <0, the scaling function | r | / <^) may prevent C (r, f, 0, t) from diverging in small distances. 
In order to illustrate Eq. (10), the Ornstein-Zernike type correlation functions Eq. (9) can be matched against it with 

Ci(r,r,0,f) = ai|r|i#i with ai = and^— (11a) 

V 9 / 

C 3 (r,f, 0 ,t) = with 03 = ^ and'^ 3 (x) = , ( 11 b) 

and = ^JvJe. All quantities are determined up to a -independent pre-factor, as one demands that all ^ -dependence 
is contained in the scaling function In both cases 77 = 0, as expected for the null-models studied. A non-vanishing 
exponent 77 is a clear signal for non-trivial long-range behavior, (i.e., when correlations on the large scale carry the 
signature of the interaction) which can therefore be considered as shaping the large scale. However, the inverse is not 
true as 77 = 0 does not necessarily mean triviality (as found in the response function for the Manna Model, Pruessner 
2013), as other correlation functions and other observables might still carry the signal of an effective long-range 
interaction even when the response function does not. The exponent 77 is normally positive, {i.e., interaction) and 
therefore fluctuations make correlations decay quicker. Beyond Tj =2 the correlations decay so quickly that coarse 
grained local degrees of freedom display Gaussian correlations (Section 2.3.4 and Pruessner (2012)). In almost all 
traditional models of equilibrium phase transitions, 77 is a small, positive quantity, with 77 = 1/4 in the 2D-Ising 
Model (Stanley 1971) being the large exception (e.g., Berges et al. 2002). 


2.1.5. Surface growth 

As an example of the use of correlation functions in the numerical detection of SOC, it is instructive to apply them 
to the study of growth phenomena closely related to SOC, such as the Edwards-Wilkinson equation mentioned above 
(Barabasi and Stanley 1995). Traditionally, exponents in the two areas have been named differently. The roughness of 
an interface 0 (r,7) above a cZ-dimensional substrate of volume V = Lf^ and linear extent L is 

= ^yd"'od"'r2 ((^(ri,7)-^(r2,f))^) . (12) 

Provided (^(ri,7)) = 0 and assuming translational invariance, this is 

w^iL,t)=Ci0,t,0,t)-^ y'dSd‘'r2C(ri-r2,7,0,7) , (13) 

an example of a sum-rule. According to Eamily and Vicsek (1985) the roughness is expected to scale like 

w^{L,t)=aL^ ’ 


(14) 
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with metric factors a and b, roughness exponent a, dynamical exponent z and universal scaling function, It is 
natural to trace the scaling of the roughness to that of the correlation function, 

C(r,t,0,f) =fl|rp“S# , (15) 

even when a number of caveats apply (Lopez 1999) in particular in the presence of boundaries or generally in finite 
systems (Pruessner 2004a). The language of interface dynamics has a long-standing tradition in SOC and a number of 
deep-running links between SOC and well understood models of surface growth have been established (Paczuski and 
Boettcher 1996; Pruessner 2003, 2012). 


Comparing Eq. (15) to the generalized form of an Ornstein-Zernike correlation function Eq. (10) implies a = 
{2 — d ~ J7)/2, which for rj = 0 reproduces the known results for the Edwards-Wilkinson equation (Krug 1997). 
Correspondingly, the correlation length is set by the growth time ^ f Eq. (15) remains valid up to a time scale set 
by the system size, t <^U-. After that, the correlation length is curbed by the system size, i.e., in Eq. (15) is replaced 
by (cf */^)). The same exponents characterizing Eq. (15) are expected to govern the two-time, two- 

point correlation function at stationarity, 

C(r,/,0,0)=fl|rr(‘'-2+r,)^|^^j ^ 

in an extension of Eq. (10). An equivalent relation is expected to hold for the response function Eq. (5). 


In the presence of a cutoff, set by the system size or other limitations, the decay of correlations on the large 
scale is characterized by the scaling function, whose typical form is that of an exponential, i.e., in Eq. (15) and 
in Eq. (16) are essentially exponentials. It is common practice to fit C(r,f,0,t) against A|r|^ exp (—r/i®) with some 
amplitude A, exponent /r and correlation length ^. The latter can be extracted very elegantly, up to the amplitude, 
by noticing that for 77 = 0 in Eq. (16) gives 0 = to leading order in <^. On a one-dimensional lattice 

(where ^ is dimensionless) this is easily verified explicitly using Eq. (11a), as 


= exp(-l/^) ^ .2 _ 1. 

(l-exp(-l/<^))2 6 


-2 


), 


(17) 


but the same proportionality holds for higher dimensions. As mentioned above, the paradigmatic form of the correla¬ 
tion function (or the propagator) in Eourier space is 


1 


(18) 


which, for small k, converges to as expected since £rC(|r|,f,0,t) is the 0-mode of the Eourier transform. Com¬ 
plicated boundary condition either spoil the structure of Eq. (18) or require orthogonal functions different from 
exp(—ikx). As such, the time separation is exemplified via the time step and the iteration, and the slower timescale 
moves with the number of external perturbations received by the system. Although all correlation functions discussed 
so far are defined on the microscopic, fast moving time scale, SOC systems normally provide a second, slow time- 
scale, whose time moves with the number of avalanches generated. Although theoretically less relevant, correlations 
have also been studied on this coarser time scale (Sokolov et al. 2014) which can be linked back to the microscopic 
dynamics (Pickering et al. 2012; Pruessner 2012). 


2.1.6. Measuring correlation functions 

There are three main reasons why correlation functions have not received much attention in experimental, numer¬ 
ical, and observational work on SOC: they require high resolution data to start with; they can be technically difficult 
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to determine e.g., (Anderson 1971); they are notoriously noisy or prohibitively expensive in terms of computational 
effort. The reason for the latter point is not least that the correlation functions have to be determined for a range of 
different coordinates ri,fi and r 2 ,f 2 to reveal the full functional dependence on these parameters. In the presence 
of boundaries, barely any of the symmetries mentioned above can be exploited to ease the computational effort. In 
the presence of translational invariance the discrete Fourier transform on a hyper-cubic lattice gives (Newman and 
Barkema 1999) 

C(k,f,0,0=E^'‘"'C(k,t,0,f) = ^(l^'(k,0P) , (19) 

r -iV 

where N = Y,v denotes the number of sites and 0(k,f) is the Fourier transform of 0(r,t) — (^(r,f)), which in the 
presence of translational invariance equals that of just 0(r,f) except for k = 0. In numerical applications, the Fourier 
transform is available as a Fast Fourier Transform (Press et al. 1992). 

Where this is computationally too expensive approximative schemes can be employed (Holm and Janke 1993) 
determining the correlation length from 1 /C(k,f ,0,f) (k^ + 1 for a few k, Eq. (18), at least for small rj. Simi¬ 

larly, taking of the correlation function numerically can produce good estimates of the correlation length, assuming 
the generalised Ornstein-Zernike form, Eq. (10), provided t] can be assumed to be small and in particular when 
d — 2 + rj =0. Up to some prefecture,the square of the correlation length is also given by the gap of the 0-mode 
C(k = 0,t,0,t) = A direct measurement of the correlations, is often hindered by the lack of symmetry. In the 
presence of conservation, SOC systems have boundaries to dissipate the energy (or particles or whatever is entering 
the system via the driving) which means that translational invariance is broken. In that case, many of the standard 
techniques fail when they rely on a standard Eourier transform. 



(a) Activity correlations in a linear-linear plot. 



plot. 


Eig. 1.— The two-point correlation function C(L/2,t,L/2-f r,f) of the activity in the Abelian version of the one¬ 
dimensional Manna Model (Manna 1991; Dhar 1999). In the language adopted in Eq. (1), ^(r,f) is the level of 
activity (i.e., at a certain point in space r and a certain microscopic time t the level of avalanching is a Poisson process 
with unit rate times the number of pairs ready to topple) measured in the middle, r 2 = L/2, and across the lattice, 
ri = L/2-|-r. (a) shows the data on a linear scale. That they collapse nicely according to Eq. (10) can be seen in 
(b), where the scaling of the abscissa is shown to be compatible with the assumption that the correlation length scales 
linearly in the system size. 
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2.1.7. Example: The Manna Model 

Figure 1 shows data of C(L/2,f,L/2 + r,f) for the Abelian, one-dimensional Manna Model (Manna 1991; Dhar 
1999) whose correlation function can be determined comparatively easily. In the Manna Model each site is occupied by 
a non-negative number of particles. As long as any site carries more than one particle, that site redistributes two of them 
among independently chosen nearest neighbors, potentially making them exceed the threshold and thereby giving rise 
to an avalanche. While the particle number is conserved in the bulk, sites toppling along the open boundary can lose 
one or two particles by moving them outside the lattice. The Manna Model is normally started from an empty lattice, 
and driven whenever the system is quiescent by depositing particles at randomly, uniformly-chosen sites. The activity 
for this model is defined in the following as the number of pairs on a site about to be re-distributed. The activity- 
activity correlation function in Figure 1 displays a long-ranged decay, whose scaling behaviour, however, becomes 
apparent only when plotted double logarithmically. In fact, the data can be collapsed acceptably well according to 
Eq. (10) with ^ = L and d — 2 + rj « 0.658, i.e., t] « 1.658, rather large compared to, say, T] = 1/4 in the Ising 
Model. Further, the scaling of the two-point activity (i.e., activity-activity) correlation function in the Manna model 
thus differs significantly from that of the propagator G, which is known to remain classical, rj =0 (i.e., of the form 
Eq. (9)), in the stationary state (Pruessner 2013). 

Various identities exist relating exponents of the activity to exponents of the avalanches (Liibeck and Heger 2003; 
Pruessner 2012, in particular p. 340). The variance of the activity density, Apa/L^, is expected to scale like 
(Liibeck 2004), which is related to C(ri,t,r 2 ,t) by the sum rule, 

^ = ZF/d"nd^r2C(ri,f,r2,f)ocL-(‘'-2+r,) ^ ^20) 

which reproduces the well known Eisher scaling law (Stanley 1971) Vx (2 — Tj) = /. In the present case //Vx = 0.41 ± 
0.04 and therefore rj = 1.59 ±0.04 (Liibeck 2004) and 2 — T] « 0.342 measured above suggests a slight mismatch, 
which might be explained by finite size effects. 



(a) Substrate correlations in a linear-linear plot. 
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(b) Attempted collapse of substrate correlations. 


10 


Eig. 2.— Similar to Eigure 1 these plots show the correlations in the inactive particles of the Manna Model (the 
substrate) measured during periods of quiescence. In the language of Eq. (1) ^(r,f) is the number of particles resting 
on r and time f, measured between avalanches, (a) suggests some short-ranged correlations, but it also indicates no 
discernible difference of these correlations for different system sizes. This is confirmed by the failure of the attempted 
collapse in (b), where the amplitudes Ai to rescale the data along the ordinate have been chosen as to facilitate the best 
collapse. 
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In contrast. Figure 2 shows the correlations in the inactive particles in the Manna Model (i.e., particles that are 
not moving around) measured during times of quiescence when no avalanche is running. While correlations do exist 
over a small number of lattice sites, the correlation length does not change with system size. This is clearly visible in 
Figure 2(a) as the data collapses without the need of any rescaling. In fact, the attempted collapse in Figure 2(b) is very 
poor and does in fact show no sign of scaling. This finding agrees with recent field-theoretical work (Pruessner 2015) 
which suggests that correlations in the substrate (i.e., the background of inactive particles) are either irrelevant or enter 
only in a very subtle way that is insignificant at large temporal and spatial scales. In other words, the substrate is an 
unsuitable place to look for correlations and SOC takes place during avalanching, not during quiescence. However, 
this finding disagrees with the traditional view that the SOC state is one of subtle correlations stored in the substrate 
(e.g., Christensen and Olami 1992; Lise 2002). Finally, we note that correlations in the substrate are mostly anti¬ 
correlations, i.e., fluctuations above the mean are repelling each other. In other words, wherever unusually many 
particles are found at one point, the environment is depleted, suggesting that the dynamics has led to a pile-up. Again, 
that ties in well with the self-organization maintaining a particular density of particles, with fluctuations only due to 
some local re-shuffling. 


2.2. Structure Functions 

The structure function provides another widely-used, two-points, statistical moment of a random variable in a 
critical system that can be used to study scaling behavior and inter-scale connections. A phenomenological analogy 
with the autocorrelation function shows that the product of field values in two points in the autocorrelation function 
is replaced by the absolute value of the increment in the definition of the structure function. The replacement offers 
an opportunity to consider various powers, q, of the increment, and thus to explore the high-order statistical moments, 
which, in turn, uncover the multifractality and intermittency properties of a system under study. The structure functions 
were first introduced by Kolmogorov (1941) (hereafter K41) in developing his turbulence theory. Note that the solar 
photospheric plasma - the medium to which a bulk of our further discussion is applied - is in a state of highly developed 
turbulence. Structure functions are defined as statistical moments of the increments of a turbulent field u(x) as 

= (|u(x + r)-u(x)|'?), (21) 

where r is a separation vector, and ^ is a real number. In the original K41 theory, u(x) is assumed to be a fluctuating 
velocity field, however the structure functions technique is applicable for any random variable, in both temporal and 
spatial domains, (e.g., Stolovitzky and Sreenivasan 1992; Consolini, et al. 1999; Buchlin et al. 2006; Uritsky et al. 
2007). For example, in Figures 3-7 the structure function technique is applied for the longitudinal component of the 
photospheric magnetic field. Structure functions, calculated within the inertial range of scales, r,(rj <r <L, where 77 
is a spatial scale where the influence of viscosity becomes significant and L is a scaling factor for the whole system) 
are described by a power law (Kolmogorov 1941; Monin and Yaglom 1975; Frisch 1995), 

5,(r)~(e,(x).r)^/3~(r)?W. (22) 

where er(x) is the energy dissipation, averaged over a sphere of size r. 

The function 1 ^ [q) describes one of the most important characteristics of a turbulent field. In order to estimate this 
function, Kolmogorov assumed that for fully developed turbulence (i.e., turbulence at high Reynolds number when the 
inertial force vastly exceeds the viscous force), the probability distribution laws of velocity increments depend only 
on the first moment, e, of the function er(x). Replacing er(x) in equation (22) by e we have 

S,{r)^{e-rfl^=C-dil\ 


(23) 
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where C is a constant. As a result, function (q) is defined as a straight line with a slope of 1 /3 

C{q)=q/3.. (24) 

Kolmogorov further realized (see also formulation of Landau’s objection concerning the original K41 theory in Frisch 
1995) that such an assumption is very rigid and turbulent state is not homogeneous across spatial scales. There is a 
greater spatial concentration of turbulent activity at smaller scales than at larger scales. This indicates that the energy 
flow and dissipation do not occur everywhere, and that the energy dissipation held should be highly inhomogeneous, 

AR NOAA 0501, Nov 19, 2003, 10:04 UT 



0.1 1 10 100 1 2 3 4 5 6 



r, Mm q 


Fig. 3.— Structure functions Sq{r) {upper left) calculated from a magnetogram of active region NOAA AR 10501 by 
Equation (21). Lower left: - hatness function F’(r) calculated from the structure functions by Equation (31) where Sq{r) 
use the longitudinal component of the magnetic held for u. Vertical dotted lines mark the interval of multifractality, 
Ar,where hatness grows as power law when r decreases. The interval Ar is also marked in upper left frame. The power 
index k is determined within Ar. The slope of Sq{r), dehned for each q within Ar, is l^{q) function (upper right), 
which is a concave for a multifractal and straight line for a monofractal. Lower right: - function h{q) is a derivative 
of {q). The interval between the maximum and minimum values of h{q) is dehned as a degree of multifractality. Ah. 
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intermittent, and follows a power law, 

where p is a real number. Then equation (22) may be rewritten as 

Sq{r) ~ (er(x) • = (er(x))'?/^ • 


(25) 


(26) 


or 

i;{q) = x{q/3)+q/3. (27) 

Equation (27) is referred to as the refined Kolmogorov’s theory of fully developed turbulence (Kolmogorov 1962a,b; 
Monin and Yaglom 1975; Frisch 1995). One can see from equation (27) that the function l^{q) deviates from the 
straight ^/3 line - the deviation is caused by the scaling properties of a field of energy dissipation. 

Important information on a turbulent field can be derived from the functions C (?) that can be obtained from 
experimental data. The value of the function at ^ = 6 deserves special attention because it defines a power index 


^^1-C(6) 


(28) 


of a spectrum of energy dissipation e(x): 




(29) 


where A: is a wave number as discussed in Section 3.1.3 below. By measuring ^(q) from experimental data and using 
equation (27) one can calculate the scaling exponent T{q/3) in equation (25) for the energy dissipation field. The 
derivative of ^(q), 


h{q) = 


dC(?) 

dq 


(30) 


can also be obtained by using the l^{q) function (Figure 3, right bottom). The deviation of h{q) from a constant value 
is a direct manifestation of intermittency in a turbulence field, which is equivalent to the term multifractality in fractal 
terminology (see further discussion in Section 3.1.3). 


2.2.1. The flatness function as an output of two structure functions 

The weakest point in the above technique is to determine the scale range, Ar, where the slope ^ (q) is to be 
calculated (see Figure 3). Abramenko (2005a) used the flatness function, defined as a ratio of the fourth statistical 
moment to the square of the second statistical moment, to visualize the range of multifractality, Ar. Another option is 
to use higher statistical moments to calculate the (hyper-)flatness, namely, the ratio of the sixth moment to the cube of 
the second: 

F{r) = S6ir)/iS2{r)Y. (31) 

For monofractal structures, the flatness, F{r) is not dependent on the scale, r. On the contrary, for a multifractal 
structure, the flatness grows as a power law, when the scale r decreases: F{r) ~ The interval Ar of the power 
law is well defined between the two cutoffs of the spectrum (see Figure 3, bottom left). The power index of the 
flatness function, K, can be used as a measure of multifractality - more complex structures have steeper F{r) spectra. 
Moreover, the interval Ar outlines the range of scales where the property of multifractality and intermittency is met. 
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2.2.2. Connection to the multifractality spectrum, f{(x) 

The function {q) is a straight line for a monofractal, due to a global scale-invariance, whereas it has a concave 
shape in case of a multifractal. The degree of concavity is usually measured by function h{q) = ^{q)/dq. All values 
of h within some range are permitted for a multifractal. For each value of h there is a monofractal with an /z-dependent 
dimension D{h) at which the scaling holds with exponent h. This representation of multifractality is based on the 
increments of the field and has its roots in the K41 theory of turbulence. A second representation is based on the 
dissipation, e, of the field energy, which relies on the K41 result stating that field increments over a distance r scale 
as (er)'/^, known as the refined similarity hypothesis Monin and Yaglom (1975). In multifractal terminology, the 
refined scaling hypothesis means that for any singularity of exponent a of er, there exists an associated singularity of 
exponent h — a/3 for the field of the same set, which has the same dimension D{h). Usually, it is very difficult to 
measure the local dissipation in the 3D space, and so one-dimensional space averages of the dissipation are typically 
used. The corresponding dimension f{a) = D{h) — (z/ — 1) is lowered by two units (for the space dimension d = 3) 
where one-dimensional cuts of a 3D structure are taken. In the literature f{a) is often referred as the multifractality 
spectrum (e.g., Feder 1988; Lawrence et al. 1993; Frisch 1995; Schroeder 2000; Conlon et al. 2008; McAteer et al. 
2010). The values of D{h), in turn, can be calculated as a Legendre transform of ^(q) Frisch (1995), 

D{h{q)) = inf^{d + qh{q) - ^{q)). (32) 

When ^(q) is concave, then for a given real value of q the extremum in Eq. 32 is attained at the unique value ho{q), 
and 

D{ho{q))=d + qho{q)-C{q))- (33) 



Fig. 4.— Structure functions Sq{r), flatness function F{r) and ^{q) function from a magnetogram of high-flaring 
NOAA AR 9077 {left. Ah = 0.48), and from a magnetogram of low-flaring NOAA AR 10061 {right. Ah = 0.06). The 
multifractality index K is the slope of F{r) calculated inside Ar. Other notations are the same as in Figure 3. 
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The result of the structure function method as applied to solar active region magnetograms are presented in 
Figure 4 (Abramenko et al. 2002; Abramenko 2005a,b; Abramenko and Yurchyshyn 2010) . The scaling behavior 
of the structure functions is different for each region. For the complex and flare-productive NOAA AR 9077 there 
is a well-defined range of scales, Ar = (4 — 23) Mm where flatness F{r) grows with the power index if = —1.17 as 
r decreases. Function l^{q) is concave and the corresponding Ah k, 0.5. This implies a multifractal structure of the 
magnetic held in this active region. To the contrary, the simple non-flaring NOAA AR 10061 (Figure 4, right) exhibits 
a flatness function that undulates around a horizontal line, which implies a monofractal character of the magnetic 
held. The function l^{q) is nearly a straight line with a vanishing value of Ah k, 0.05. Time profiles of Ah for the 
two active regions are compared in Figure 5. The non-flaring NOAA AR 10061 persistently displays lower degree of 
multifractality, as well as lower X-ray flux, than the flaring NOAA AR 9077 does. Figure 6 demonstrates the statistical 
relationship between the multifractality index, k and a flaring index, A for 214 regions (Abramenko and Yurchyshyn 



time, hours 


Fig. 5.— Time variations of the measure of multifractality, A/i (left axis), and GOES soft X-ray flux (right axis, dashed 
lines) plotted for six-hour time intervals for the two active regions. Data for NOAA AR 9077 (red lines) were obtained 
between 17:00 and 23:00 UT on July 13, 2000 and data for NOAA AR 10061 (green lines) refer to an interval between 
11:00 and 17:00 UT on August 9, 2002. 





-17- 


2010 ), from which it is clear that the higher degree of multifractality of the magnetic field may be associated with 
stronger flare productivity of an active region. Here the flare index A characterizes the flare productivity of an active 
region per day, being equal to 1 (100) when the specific flare productivity is one Cl.O (Xl.O) flare per day. More 
examples of multifractality spectra f{a) are shown in Figure 7 (Abramenko and Yurchyshyn 2010). One can see that 
the most complex and flare-productive regions (left frame in Figure 7) exhibit broader spectra as compared to that of 
non-flaring regions (right frame). This means that a set of monofractals that form an observed multifractal, is much 
more broad in flare-productive regions as compared to non-flaring regions. 



- 2-10 1 2 
Multifractality Index, k 


Fig. 6.— Flaring index. A, plotted versus the multifractality index, K, for 214 regions. The Pearson correlation 
coefficient is -0.63. From Abramenko and Yurchyshyn (2010).) 
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(a) high-flare productivity 


(b) low-flare productivity 


Fig. 7.— Multifractality spectra, f(a), plotted; /eft - for regions of high flare productivity; rig/it - for regions of low 
flare productivity and for a plage area. Spectra on the left frame are more broad than that on the right frame. 
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2.3. Application Oriented Methods 

As discussed in Section 2.1, the classical autocorrelation SOC detection methods are explicit in theory, but are 
often challenging in terms of practical application to physical systems, such as the solar atmosphere or the tectonic 
environment. Over the previous 25 years, and through the evolution of several numerical SOC models created to ex¬ 
plain existing physical systems, a variety of application-oriented methods have been developed that together comprise 
a useful toolkit for the detection of the SOC state. In principle, when the SOC state is reached the system experiences 
instabilities of all sizes, clustered in cascades of elementary events, or avalanches, all triggered by fixed and small 
(with respect to the critical threshold) or variable, but statistically small, perturbations (for the latter set of SOC mod¬ 
els, see Georgoulis and Vlahos (1996, 1998)). The main feature of this marginally stable SOC state, where a given 
small perturbation can cause avalanches of all sizes (e.g., Kadanoff 1991; Newman etal. 1996) is precisely the absence 
of a preferred scale for avalanche size. This leads to robust power laws if one examines the distribution function of 
the event sizes (Section 3.2.1). In this sense, a nonlinear dynamical system realizes the SOC state as a statistically- 
stationary state far from equilibrium. We review these two attributes of marginal stability and statistical stationarity 
as practical detection methods for an SOC state. We then present a recent non-evolutionary diagnostic SOC-state test 
and finally discuss block-scaling methodology is detail. 


2.3.1. Marginal stability: a spatially averaged critical quantity 

The diagnostic SOC detection method of marginal stability is based on the stabilization of a spatially averaged 
system parameter, i.e., the parameter compared with the critical threshold. Applying this method to the classical 2D 
cellular automaton sandpile model of Bak et al. (1987), it is assumed that each point i = (x,y) of the square grid 
corresponds to the space occupied by a sand grain. The field variables in this model are the height h{x,y,t) and the 
slope G{x,y,t) of the accumulated sand at every point, i = {x,y) of the system and in every time step f, of its evolution. 
Referring to the classical cellular automaton, both space and time are discretized; the automaton consists of a discrete 
grid , e.g., (x,y) in 2D, where each grid site has a position vector i with integer components. The automaton also has 
two discrete time-scales, namely an integer time step, t, that increases by one with each application of the automaton 
rules, and an integer iteration that increases by one each time the system is perturbed. The slope G{x,y,t) at a specific 
point of this automaton’s sandpile and for the specific time t is defined as the height difference between the height 
h{x,y,t) at the point i = {x,y) and the average height of the adjacent grid points h{t ), 

h{t) = -^[h{x+ l,yd)+h(x— \,yd)+h{x,y+ \,t) +h{x,y— 1,/)] . (34) 

Therefore, the slope G{x,y,t) is defined as G{x,y,t) = h{x,y,t) — h{t). The transition rules describing the evolution 
of the system when a sand grain is added at a random point i = {x,y) of the grid at time t are defined as h{x,y,t) —^ 
h{x,y,t) + 1- The instability criterion embedded in the transition rules of the system reflects a critical value of the 
slope Gc. A point i = {x,y) of the system is considered unstable when the inequality G{x,y,t) > Gc is fulfilled. When 
such an instability occurs at the point i = {x,y) and at the time f, then the dynamical system responds at the time f -f 1 
according to the following evolution or redistribution rules. 


h{x,y,t -f 1) = h{x,y,t) — 4 , 

(35) 

h{x± \,y,t+ 1 ) = h{x3i l,y,t) + 1 , 

(36) 

h{x,y± l,t-|-l) = /z(x,y±l,f)-|-l . 

(37) 


Transition and evolution rules comprise the driving and relaxation mechanisms, respectively, that inexorably lead the 
system to marginal stability. A practical SOC-state detection mechanism based on this marginal stability reached by 
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system in such a state was presented by Georgoulis (2000). This mechanism monitored the temporal evolution of 
the mean value of the field variable(s) that determine(s) the instability threshold for the system. For the Bak et al. 
(1987) model described above, Georgoulis (2000) monitored the temporal evolution of the mean height H{t) of the 
sandpile throughout the grid, where H{t) = with i being the position vector. Equivalently, one can monitor 

the temporal evolution of the mean slope G{t) throughout the grid, where G(f) = ^ ^ as SOC can be reached in 
both critical-slope and critical-height cellular automata models (Kadanoff et al. 1989). 

Figure 8 presents the temporal evolution of the mean height H{t) and the mean slope G{t) for a 3D sandpile 
cellular automaton model with dimensions 20 x 20 x 20. Initially both the mean height H{t) and the mean slope 
G(f) are increasing. This ascending course corresponds to the sequence of the metastable states, through which the 
system evolves towards the SOC state. This marginally stable state is reflected in the stabilization of both variables 
after the dashed vertical line. This line determines the time, in system iterations, after which the system enters the 
SOC state, generating avalanches lacking a characteristic scale in size or duration. Figure 8 also shows that after the 
SOC state is reached, the mean slope G(f) stabilizes around a value slightly lower than that of the critical threshold 
Gc- In the cellular automaton model used in this example, the critical threshold (horizontal dashed line) is Gc = 10, 
in arbitrary system units. In addition, the SOC state is reached after ~ 4.8 x 10® iterations, which corresponds to 



Time (iterations x 10^) 


Fig. 8.— Time evolution of the mean height H{t) and the mean slope Git) for a 3D statistical Flare cellular automaton 
sandpile with dimensions 20 x 20 x 20. The dashed vertical line corresponds to the time at which the system enters 
the SOC state. The horizontal dashed line corresponds to the critical threshold value of the slope, which defines the 
instability criterion of the system. From Georgoulis (2000)). 
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^ 0.6N^ iterations, where N = 8 x 10^ is the number of nodes, or grid sites, in the SOC system used here. This 
number of iterations is in order-of-magnitude agreement with the prediction of Charbonneau et al. (2001) regarding 
the number of iterations needed to reach SOC (~ N‘^, where d is the Euclidean dimension of the system), although the 
proportionality factor here is ~ 1, where in the prediction of Charbonneau et al. (2001) it is typically ^ 1. Possibly 
this is due to the fact that the statistical flare model of Georgoulis and Vlahos (1996, 1998), which is the one used in 
Figure 8, does not apply a fixed, infinitesimal driving, but rather uses a perturbation of variable amplitude that is small 
on average as compared to the critical threshold. This appears to shorten the driving time needed for the system to 
reach the SOC state. 

The same method was adopted by Dimitropoulou et al. (2011) for the detection of the SOC state in a 3D cellular 
automaton that included vector, rather than scalar, magnetic fields such as the seminal models of Lu and Hamilton 
(1991) and Lu et al. (1993). The novel element of this work, however, is that the magnetic field vector is data- 
driven, i.e., relying on actual solar active regions. The model uses an observed photospheric vector magnetogram of 
a given active region and extrapolates it via a nonlinear force-free extrapolation (Wiegelmann 2008) into the overlay¬ 
ing corona, thus obtaining the initial 3D vector field. The configuration is subsequently evolved into the SOC state 
using conventional cellular-automata rules. This model has been coined the static integrated flare model (S-IFM) by 
Dimitropoulou et al. (2011) because it refers to a single, simultaneous magnetogram. In this model it is assumed that 
instabilities occur if the magnetic field stress exceeds a critical threshold. For every site r within a cubic grid with 
dimensions 32 x 32 x 32, the magnetic field stress Gav(r) is calculated as Gav{r) = |Gav(r)| where 

G«v(r)=B(r)-—^B„„(r) , (38) 

«« ™ 

where nn is the number of nearest neighbors for each site r and B„„(r) is the magnetic field vector of these neighbors. 



Fig. 9.— Average Laplacian Gav over the grid for 3 x 10^ timesteps for NOAA 10570. Gav increases gradually for 
1.4 X 10^ timesteps, after which the SOC state is reached, with Gav ^ Gcr = lOG. From Dimitropoulou et al. (2011). 
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Depending on the location of each site within the volume, the number of nearest neighbors nn can be 3, 4, 5, or 6 in 
3D, for an edge, vertex, boundary or interior location of the examined grid site, respectively. As Gav is related to the 
diffusive term of the induction equation, it was selected by Dimitropoulou et al. (2011) to be compared against the 
critical quantity of the system such that every site r = {i,j,k) for which the inequality Gavjjf. > Gcr = lOG is satisfied 
is considered unstable and undergoes magnetic field restructuring according to specific evolution rules. By monitoring 
the volume average Gav of the critical quantity Gav, it was shown that Gav increases gradually during the continuous 
driving of the system. When the system reaches the SOC state, Gav stabilizes around a value slightly lower than the 
threshold value Gcr- Figure 9 shows Gav value over 3 x 10^ time steps for a solar active region (NOAA AR 10570). 
Gav increases up to time step ~ 1.4 x 10^, thereafter asymptotically tending to the critical threshold at Gcr = lOG. 
A second indication that the system has reached the SOC state is that the total volume energy attains an asymptotic 
value stemming from the competing tendencies of injecting energy in the system via driving and dissipating it via 
relaxation events. Figure 10 shows the logarithm of the volume magnetic energy Etotaft after each scan of the grid for 
possible re-distributions. Etotaft shows when the system appears to reach the SOC state, namely at ~ 10® iterations, or 
~ (1 /32) X where N = 32^ is the number of system nodes in this case. This is again dimensionally consistent with 
the prediction of Charbonneau et al. (2001), although the proportionality factor is much smaller than the one predicted 
in that study, even though the driving perturbations in Dimitropoulou et al. (2011) have a fixed amplitude. 



Fig. 10.— Total Volume Energy, \ogi(i{Etotaft), after each redistribution for NOAA AR 10570. As in Figure 9, Etotaft 
increases gradually until an asymptotic stable state is reached. From Dimitropoulou et al. (2011). 
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2.3.2. Statistical stationarity: number of avalanches per fixed time interval 

Statistical stationarity can also be used as an applied diagnostic method towards the detection of the SOC state. 
This is based on the premise that after a dynamical system has entered the SOC state, the number of avalanches 
produced within a fixed time interval will vary around a well defined average value Georgoulis (2000). Figure 11 shows 
an example of this variance that corresponds to the same 3D cellular automaton model of Georgoulis (2000) described 
in Figure 8. In particular, Figure 11 shows a time series of the number of avalanches produced in fixed time intervals 
consisting of 1000 model iterations. A new iteration is triggered when a sand grain is added to the modeled sandpile at 
one specific, randomly chosen, grid point {i.e., h{x,y,t) h{x,y,t) + 1, as above). In accordance to conventional SOC 
models, the driving of the system is not continuous, with each new iteration requiring the complete relaxation of all 
avalanches in the system. As a result of the statistical stationarity embedded in the SOC state dynamics, the number 
of avalanches per 1000 iterations varies around a well defined average value of ~50 events, regardless of event size. 

The same method was applied to the static, data-driven, integrated flare model (Dimitropoulou et al. 2011), as 
described in the previous paragraph. Figure 12 shows the average number of avalanches, this time for a single vector 
magnetogram of the observed NOAA AR 11158, as a function of the simulation iterations. The driving of the system 
is also not continuous and is applied to a single, random grid point as long as there are no ongoing avalanches. It is 
shown that after approximately the first 130,000 iterations the average number of the produced avalanches stabilizes 
around ^450 events per 1000 iterations, which attests to the statistically stationary SOC state reached by the system. 


2.3.3. Non-evolutionary diagnostic SOC-state test 

A third SOC-state test is made possible from the coupling between two data-driven solar flare cellular automata 
models: the static (S-IFM) model and the dynamic (D-IFM) model. Rather than detecting the SOC state in line with 
the previous tests (i.e., on an evolution time series of a possible SOC system), this non-evolutionary diagnostic aims 
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Fig. 11.— Time series of the number of avalanches produced per 1000 iterations for the same 3D statistical flare 
cellular automaton model discussed in Figure 8. A statistical stabilization of the average number of events is shown, 
after the system has reached the SOC state, beyond the first 2 x 10® iterations. From Georgoulis (2000). 
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to determine whether a given 3D snapshot magnetic configuration could be in the SOC state. Both the classical (e.g., 
autocorrelation test of Section 2.1) and the applied methods of marginal stability and statistical stationarity tests rely 
on an SOC-state detection based on a continuous monitoring of the evolution of a potential SOC system. This non¬ 
evolutionary test instead offers an indication of whether an instantaneously observed system is possibly in an SOC 
state, among other possible physical mechanisms that may have led it to the observed configuration. 

A brief description of the D-IFM method is attempted here for context: in D-IFM, the single vector magnetogram 
of S-IFM is replaced by a time series of vector magnetograms of a given active region. Each magnetogram of the time 
series is subjected to the S-IFM methodology, i.e., an initial nonlinear force-free extrapolation to obtain the 3D coronal 
magnetic field and a randomly driven evolution into the SOC state. Each magnetic configuration is confirmed to have 
reached the SOC state through the marginal stability and statistical stationarity tests. The D-IEM then proceeds by 
slowly driving the magnetic configuration from the one 3D SOC snapshot to the next via a spline interpolation of the 
magnetic field components. The number of iterations is typically >> 1 for observational cadence of the order tens of 
minutes and depends on the Alfven time required to cross a distance equal to the line element (pixel size) assuming a 
constant, typical coronal Alfven speed of 10^ cm/s (Dimitropoulou et al. 2013, Table 3). In this course, avalanches 
occur and are relaxed, giving rise to a sequence of SOC-state events with properties that are studied statistically. 
Eigure 13 depicts this basic D-IEM concept applied to a time series of 7 vector magnetograms of the observed NOAA 
AR 8210. Avalanches occur when the critical threshold of the magnetic field Laplacian is exceeded. Moreover, 
numerous sequences, or groups, of 3D configurations can be obtained, for each of which one may independently apply 
the D-IEM and collect the statistics jointly. 

It is this coupling between the static and dynamic models that inspires the concept of the following non-evolutionary 
diagnostic SOC-test. The principal idea is to apply the S-IFM to an observation (vector magnetogram), leading the 
initial NEFF field solution into a SOC-state magnetic configuration. The random forcing of the S-IFM will give rise 
to a very different SOC-state configuration, as compared to the initial NEFF field solution. Then, the same instability 
criterion is used to revert the configuration to the initial NEFF field solution via the D-IFM, i.e., through a continuous 



iterations 

Fig. 12.— Time series of the average number of avalanches produced per 1000 iterations for the static, data-driven 
cellular automaton of NOAA AR 11158. A statistical stabilization of the average number of events is shown, after the 
system has reached the SOC state, beyond the first 130,000 iterations. From Dimitropoulou et al. (2011). 
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interpolation. Since the final S-IFM snapshot is proved to have reached the SOC state and the D-IFM demonstrably 
retains the SOC characteristics, reverting this snapshot to the original NLFF field solution via the D-IFM is a good 
indication that the initial NLFF field is indeed in a SOC state. This would be impossible to claim otherwise for any 
given static 3D magnetic field solution. 

Figure 14 presents the S-IFM part of the non-evolutionary diagnostic SOC-test concept applied to the observed 
NOAA AR 11158. Figure 14a depicts the vertical component of the studied photospheric vector magnetogram, while 
Figure 14b shows the preprocessing necessary in order to apply the S-IFM, namely the re-binning of the magnetogram 
into a grid of 32x32 (left) and the subsequent 3D nonlinear force-free extrapolation (right). The choice of coarse grid 
resolution is determined by computational power available for the iterations. Figure 14c illustrates the photospheric 
vertical field component (left) and the corresponding 3D coronal configuration (right) after the S-IFM application 
for 2.5 X 10^ iterations. Notice the severe distortion of the magnetic field vector, caused by the randomness of the 
S-IFM forcing. This configuration, however, is both a valid (i.e., divergence-free) magnetic field solution and is 
demonstrably in the SOC state. Retaining the same instability threshold, the D-IFM is then applied, aiming to revert 
the 3D configuration of Figure 14c into that of Figure 14b, with the results shown in Figure 15. Evidently, the system 
reverts back to the configuration of Figure 14b after ~ 10^ iterations. The marginal stability test shows that the system 
remains in the SOC state until the end of the simulation, and therefore in the course of the continuous interpolation 



Fig. 13.— Graphical description of the D-IFM, applied to 7 vector magnetograms of NOAA AR 8210: each vertical 
sequence indicates a separate application of the S-IFM to a single IVM vector magnetogram. This leads to 7 3D SOC- 
state magnetic field configurations that can be evolved indefinitely. For each horizontal group of 7 3D configurations, a 
spline interpolation progresses the magnetic field vector from the one to the next configuration, collecting avalanches 
and their properties. This action corresponds to a single application of the D-IFM. In this test, 16,235 events (i.e., 
septuplet groups) have been collected in order to attain sufficient statistics. From Dimitropoulou et al. (2013). 
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Fig. 14.— a) Observed vertical component of NOAA AR 11158 (13 February 2011, 15:58:12 UT). b) Left: re¬ 
binned photospheric magnetic field to grid dimensions 32 x 32. Right: extrapolated coronal magnetic field with grid 
dimensions 32 x 32 x 32. c) Similar to b, but after the S-IFM application for 2.5 x 10^ iterations. 



Fig. 15.— Non-evolutionary SOC test run on a snapshot of the observed NOAA AR 11158, shown in Figure 14. The 
S-IFM has brought the snapshot to an SOC state after ^ 0.8 x 10^ iterations. This is confirmed by the stabilization of 
the averaged slope in the grid (curve). To ensure the unambiguous evolution to the SOC state, the S-IFM is applied 
for an additional 2.5 x 10^ iterations (vertical line). The system is then reverted back to the initial 3D-extrapolated 
magnetic configuration via the D-IFM, reaching it after ^10^ iterations, without exiting the SOC state. 
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to the initial 3D field. Continuous interpolation would not be possible if the critical threshold Gcr for the extrapolated 
field in Figure 14 was not the same with the one in the S-IFM. 

This simple SOC diagnostic suggests that both observed and force-free extrapolated solar magnetic configura¬ 
tions may already be in an SOC state - at least this is indicated by the successful test on NOAA AR 11158. It should 
be followed by including an investigation on how a far-from-equilibrium, SOC, state prevails on a force-free equilib¬ 
rium magnetic field solution. If confirmed this finding may have important ramifications on whether the global solar 
magnetic field, at least the low-j3 corona, is into an SOC state or whether this feature restricts to (many, most, or all) 
active regions. This aligns with the discussion on open problems and questions in SOC applicability, detailed in the 
review of Aschwanden et al. (2014). 


2.3.4. Block scaling 


A sum rule, similar to equation (20) above relating Apa and C (r i , f , r 2 , f ), can be used to extract scaling in systems 
when very little data is available. Although the basic concept also applies to the variance and thus to the two-point 
correlation function, it can be applied much more directly to one-point functions, i.e., to the basic degree of freedom 
0(r,f) (the local activity, energy density, particle density etc.). SOC occurs only right at the critical point, therefore 
the globally averaged activity (the order parameter) is normally very small. Although there are strong spatio-temporal 
fluctuations {i.e., the activity might flare up locally and even globally on occasion) the local activity (or generally order 
parameter density) can be averaged spatially over local patches. In the following section, these patches are referred to 
as blocks. There are N = [Ljtf such blocks of linear extension f in a cZ-dimensional system, V, with overall linear 
extension, L. Within each such block B, (such as illustrated in Figure 16a) the local activity density can be defined as 

= (39) 

first suggested by Binder (1981) for the order parameter in a ferromagnetic phase transition. Obviously, the arithmetic 
mean over the blocks is invariant under a change of £, because 

= ^ ^ ( 40 ) 

independently of £. One may introduce, however, a level of activity T, effectively a threshold, which has to be present 
somewhere in the patch if the patch is to be considered active, say 


= 0(max{0(r,t)|r G B,} — T) , 


(41) 


where 9 denotes the Heaviside theta function and max{0(r,f)|r G B,}) is the maximum activity 0(r,t) in the block B,. 
As a result is unity if ^(r,f) exceeds T somewhere in an active block. Otherwise, it vanishes. To facilitate better 
data analysis, 0(r,t) may be a function of the original raw data, with a background subtracted and/or the modulus 
taken to make it non-negative. Conditioning the average to active blocks produces the conditional activity 






(42) 


i.e., p is the average activity exceeding the threshold. This quantity displays a dependence on £, as opposed to Eq. (40) 
(which displays so such dependence). In the presence of correlations, non-vanishing a, (f) is indicative of large levels 
of activity in the whole block, such that p{t,£) should increase as £ decreases. This is strictly true for T = 0 and 
non-negative ^(r,t), in which case = 'Tii because = 0 implies 0,(f) = 0 if T = 0. In that case 
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a, (f)/A^ cannot increase as £ decreases and so p{t,£) increases with decreasing it is a matter of standard finite size 
scaling that p{t,£) £-Plv± (Pruessner 2008) with j3/Vx = (t/ — 2 + i])/2 from the usual scaling relations (Liibeck 
2004; Pruessner 2012). \f T = 0, the scaling is driven by the dominator in Eq. (42) and amounts to counting the 
number of blocks containing a certain level of activity ^ (r, f). The procedure is then not dissimilar to the box-counting 
method used in the study of fractals (Falconer 2003; McAteer et al. 2005). 

The same behavior, p {t,£) expected for T > 0, as the fraction of blocks with some activity above 

the threshold decreases with decreasing £, while those blocks i containing such high levels generally have a higher 
average activity Pi{t), i.e., the numerator is expected to increase and the denominator to decrease with decreasing i. 
As suggested by the exponent 77 , see Eq. (10), the scaling of p{t,£) is indicative of correlations. If blocks are large, 
then most of them will exceed the threshold somewhere (i.e., they will be active) and in fact p{t,£) approaches the 
unconditional average Eq. (40) as f — 7 L (as long as the threshold is smaller than the global maximum). If ^(r,f) were 
completely independent at different r then selecting them according to activity exceeding a threshold amounts to a 
random, independent selection. Provided only that the blocks are big enough that the single site where the threshold 
is exceeded does not introduce a significant bias, p{t,£) will barely increase with decreasing £, even when working on 
a lattice. Correlations, however, have the effect that regions with an activity beyond a certain threshold are generally 
more active, or, in the case of anti-correlations, significantly less active. 

A relation similar to p(t,£) the variance of the conditional activity, that is the variance 

of 0,(f) conditional to 0(r,f) exceeding some threshold within the block. In effect, block scaling gives access to finite 
size scaling, without changing the system size. In block scaling, the cutoff in correlations, avalanche size distributions 
etc., is implemented not by the system size, but by the block size. However, the linear extent of the block £ is an 
additional scale whose upper cutoff is set by the system size. Proper asymptotic scaling can be expected only when 
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(b) Block scaling of the conditional activity extracted 
from the image. 

Fig. 16.— A block scaling analysis of a snapshot of an HMI Magnetogram (11 Feb 2014). (a) The large quadratic patch 
covering most of the sun (2560 x 2560 pixels) is divided into smaller blocks (here 5x5 blocks of linear extension 512, 
some of which are labelled) as shown by the dotted lines, i.e., L = 2560, f = 512. (b) Processing the data as described 
in the text produces a narrow scaling region with an approximate exponent 0.66. Ordinate p{t,£) denotes the activity 
at the time t when the snapshot in (a) was taken averaged over those blocks of size £ which exceed a (high) threshold 
T somewhere within the block. 




(a) HMI Magnetogram from 11 Feb 2014. 
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IjL^ 1. On the other hand, a (the lattice spacing or some other microscopic cutoff) must be fulfilled to avoid 
some smaller scale physics or other effects such as resolution limitations to take over and dominate the behavior of 
p Block scaling therefore is a form of intermediate scaling (Barenblatt 1996). 

Nevertheless, the block scaling method provides access to a whole range of scales, even when, ultimately, it can¬ 
not replace finite size scaling. It is a tool to quantify correlations allowing a possible universality class to be identified. 
It has the advantage of requiring little data such as a single but highly resolved snapshot. It is in effect a sub-sampling 
scheme (Havlin and Bunde 1996), designed to extract as much information as possible from a (comparatively) sparse 
source. However, although block scaling instantly indicates the presence of correlations and its scaling, it cannot serve 
as an unique indicator for the presence of SOC. 

Figure 16 shows the results of a block scaling procedure applied to a full disk solar magnetogram. By design, 
this process specifically filters the active region patches from the quiet Sun. The data encoded in the grey-level of 
the magnetogram were processed by taking the modulus of the deviation from the overall average, and considering as 
active only those regions which are close to the maximum. In other words, the magnetic field in regions that count 
as active deviate very strongly from the mean magnetic field. Figure 16b shows a narrow region of power law, which 
may terminate or bend for very small patch sizes, where the analysis gets close to the resolution limit. Correlations of 
strong active regions are of course expected and Figure 16b shows j3/v « 0.66 and therefore 77 « 1.32 in the present 
case, again comparatively large. For comparison, 77 « 1.54 in the Manna Model (Liibeck 2004; Pruessner 2012) in 
two dimensions. 
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3. Detection of SOC-state events 

With a powerful set of tools designed to study the correlations expected to be present between features in SOC 
systems, we now turn our focus to the question of what determines a feature. In this context a feature is considered 
as collection of density enhancements in space, a variation in time, or a variation of density enhancements in spatio- 
temporal data. In this section we discuss the relevant problems with each method, and review some method-specific 
tools that have been determined as useful tools for analyzing SOC systems. 


3.1. Feature Detection in the Spatial Domain 

3.1.1. Thresholding 

Feature detection in space usually consists of dealing with a 2-dimensional greyscale image captured on a charge- 
coupled device (CCD), and often calibrated {e.g., simple CCD considerations of flatfielding, dark subtracting, etc. have 
been removed). However, these data still remain in digital number (DN) space. As such, the scientist usually considers 
a series of image processing routines, {e.g., based on standard procedures available in Falconer and Woods (2008), or 
Starck and Murtagh (2006)) that can be used to identify potential SOC features, to separate them from any noise or 
non-SOC background, and to characterize them for further analysis. One of the simplest approaches is to apply a fixed 
threshold in DN space, and group contiguous pixels into one feature. One of the earliest uses of this thresholding and 
grouping was in studies of colloidal dynamics or Brownian motion (Perrin 1920; Crocker and Grier 1996), and the 
use of such an algorithm extends to diffusion limited aggregation (Efron 1982), particles in Saturn’s rings (Zebker 
et al. 1985), and urban growth (Batty et al. 1989). The case study of solar bright points - small scale, short lived 
brightenings in the solar corona - provide some insight into the power of such a method. The threshold is usually 
considered at 2 or 3 standard deviation amplitudes above a background mean {e.g., McAteer et al. 2002, 2003). By 
adding on rules regarding feature size and feature lifetime (McAteer 2003), this procedure makes it possible to track 
features over a sequence of images {e.g., DeForest et al. 2007; Lamb et al. 2008, 2010; Kirk et al. 2012, 2013). With 
such set of extracted features, the final step is a search for correlations and power laws in their distributions (Krucker 
and Benz 1998; Parnell and Jupp 2000; Parnell et al. 2009). Although thresholding and grouping provides a simple 
and convenient method of identifying features, it is also prone to problems with sensitivity in the chosen threshold and 
in differentiating between feature disappearance and feature clumping. 


3.1.2. A volumetric consideration 

One method to overcome the known problems associated with thresholding and grouping is to use multiple im¬ 
ages of the same feature, as observed at different wavelengths. In astrophysical observations, power-law distributions 
of fluxes or fluences of candidate SOC events have been measured in almost every wavelength, from gamma-rays, 
hard X-rays, soft X-rays, EUV, visible light, to radio wavelengths. While numerical lattice simulations of SOC mod¬ 
els quantify the size of an SOC event simply by the number of active nodes that are unstable and subject to a local 
re-distribution during any time of an SOC avalanche, the size of an astrophysical SOC avalanche can only be quanti¬ 
fied in terms of an observed flux or fluence i.e., the time-integrated flux over the duration of an avalanche. However, 
astrophysical fluxes or intensities, with physical units of energy per time unit, are wavelength-dependent, and thus 
depend on the instrumental wavelength filter response function, expressed as a function of emission measure per tem¬ 
perature unit, R{T). There are different methods to convert the observed flux into wavelength-independent quantities 
that can be suitable for the characterization of the size of an SOC avalanche: conversion into radiated energy, i.e.. 
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E = riphathy = ripi^othclX, where Uphot is the number of photons that produce a flux Fx\ conversion into an emission 
measure by inversion of the flux Fx = j[dEM/dT] R{T)dT ; conversion into thermal energy = ^neksTeV, which 
requires a determination of the electron density {e.g., from the volumetric emission measure, He = y/EM/V) and the 
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Fig. 17.— Seven multi-wavelength EUV images of the X2.2-class flare observed with AIA/SDO on 2011-Feb-15 
01:50:00 UT, in the wavelengths of94A, 13lA, 17lA, 193A, 21lA, 304A, and 335A. The spatial scale of an image 
side is 0.3 solar radius 200 Mm) and the flare area is indicated with black contours at the 50% and 75% peak 
flux level. 
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electron temperature Tg. Whatever quantity is preferred to characterize the size of an SOC avalanche, this is an extra 
step that is usually not part of any numerical or mathematical SOC theory. 



Fig. 18.— Size distribution of 155 solar flare areas, obtained in 7 different wavelengths (represented in different colors 
with the wavelengths indicated on the left side, and the index of the power law, Ua, on the right side). The 155 flare 
events include all M- and X-class flares observed with AIA/SDO during 2010 May 13 and 2012 March 31. From 
Aschwanden ef al. (2013). 

A study of how events appear in different wavelengths provides insight on the spatial structuring of an SOC 
system. Figure 17 shows 7 EUV wavelength images of a large solar flare just at the peak of the emission, observed 
with AIA/SDO on 2011 February 15, 01:50 UT. A bright sigmoidal white structure is evident in the core of the active 
region, evidence of a high emission measure and a high-density heated plasma, confined in a helically twisted magnetic 
filament. Brightness contour levels at 50% and 75% of the flux maximum, include somewhat less dense heated plasma 
loops that sutTound the core, and make up a substantial fraction of the active region. Using 50% contours to demarcate 
the flare area A(f), the relative size varies considerably across different wavelengths, with a minimum size in the 94 
A Alter, and a maximum size in the 193 A Alter. To measure the actual flare area A (f), one has to subtract a pre-event 
background image A(fo), which will Alter out all static emission from the active region. It is usually not possible to 
know a priori which wavelength is the best to measure the flare area, or what flux threshold level is most appropriate 
to define the flare area. Thus, it is advisable to measure the flare area with different thresholds and in different 
wavelengths, in order to determine any possible nonlinear scaling between different wavelengths, which could in turn 
affect the slope of the power-law distributions of flare areas, N{A). Such a study has been performed with 5 different 
threshold levels and 7 wavelength Alters for 155 flares (Aschwanden et al. 2013). The resulting flare area distributions 
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are shown in Figure 18, after normalizing the flare area to the same flux threshold. Almost identical indices are obtained 
for the flare areas obtained in the 7 wavelength filters in Figure 18, which indicates that the flare areas measured in 
different wavelengths are statistically either identical or differ only by a fixed proportionality constant. The individual 
indices are also tabulated in Table 1. This result simplifies future analysis enormously, because it essentially implies 
that the choice of wavelength does not affect the statistical distributions of geometric parameters, such as the size 
distribution of lengths L, areas A, or volumes V of candidate SOC events. 



GOES flux [W/m^] GOES flux [W/m"] GOES flux [W/m^] GOES flux [W/m"] 


Fig. 19.— Correlations between the observed fluxes in 7 different AIA wavelengths with the GOES flux Fgoes for 
155 M- and X-class flares observed with AIA/SDO. From Aschwanden and Shimizu (2013). 

A complementary study of the wavelength dependence of observed fluxes provides further insight into SOC pro¬ 
cesses. Figure 19 shows scatterplots of the 7 AIA flare peak EUV fluxes with the higher energy (GOES) soft X-ray 
flux, for the same set of 155 M- and X-class flares (Aschwanden and Shimizu 2013). Apparently there exists a correla¬ 
tion between each of the EUV fluxes and the soft X-ray flux. The cross-correlation coefficients vary from CCC = 0.82 
for the 193 A filter, which shows the closest correlation with the GOES 1-8 A flux due to their overlapping high- 
temperature response (i.e., the 193 A filter is sensitive to the Ee XXV line at a temperature of Tg m 20 MK), down to 
CCC — 0.48 for the 304 A filter, which is most sensitive to cooler chromospheric plasma. Although the proportion¬ 
alities between the EUV and soft X-ray fluxes have some significant scatter, their size distributions are similar, as the 
indices ap listed in Table 1 demonstrate. Consequently, it is reasonable to also expect near-proportionality for linear 
regression fits between the EUV and soft X-ray fluxes, i.e., Fpijv ^sxw ^ scaling exponent of 7 « 1. Indeed, Ta¬ 
ble 1 shows an average exponent of 7 = 1.1 ± 0.2 for these 7 wavelengths. This important result of near-proportionality 
of EUV to SXR fluxes implies the wavelength independence of flux size distributions, which again eases comparisons 
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Table 1: The index of the power law of size distributions of flare areas a a and fluxes Fx, and scaling exponents 7 
(discussed in Figure 19) for 155 flares observed with AIA/SDO observed in 7 wavelengths. 


Instrument 

Wavelength 

A [A] 

Index of 
power law of 

area 

CCa 

Index of 
power law 
flux 

ap 

Cross correlation 

exponent 

AIA vs. GOES 

7 

AIA 

94 

2 . 0 ± 0.1 

2.2±0.04 

1 . 02 ± 0.12 

AIA 

131 

2 . 2 ± 0.2 

2 . 0 ± 0.02 

1 . 10 ± 0.12 

AIA 

171 

2.1±0.5 

2 . 0 ± 0.1 

0.90±0.11 

AIA 

193 

2.0±0.3 

2 . 0 ± 0.1 

1.19±0.10 

AIA 

211 

2.0±0.4 

2 . 1 ± 0.2 

0.87±0.08 

AIA 

304 

2 . 1 ± 0.2 

2.1±0.9 

1.36±0.25 

AIA 

335 

1.9±0.2 

1.9±0.1 

1.17±0.13 

GOES 

1-8 


1.92 


FD-DOC prediction 


2.00 

2.00 

1.00 


of SOC statistics in astrophysical objects considerably. 


3.1.3. Turbulence and Fractals: A direct 2D fingerprint of 3D SOC? 

Direct imaging has the potential to provide a direct fingerprint of detecting SOC in the spatial domain. Under 
this paradigm, it is assumed that any SOC system will involve power laws across spatial scales, and that this will 
manifest in terms of turbulence and fractality (McAteer et al. 2010; McAteer 2013, 2015). Indeed, since Kolmogorov 
(1941) and Mandelbrot (1975) first introduced the ideas of turbulence and fractals, respectively, complex systems have 
been found to be ubiquitous in many areas of human and natural sciences. Spatial power laws provide the connection 
between turbulence and SOC as discussed above in Section 2.2. The calculation of the spatial energy spectrum is given 
as 

E{k)^k-P, (43) 

where the spatial energy, E, varies with wavenumber, k, risen to a scaling index, j3. (where j3 = 5/3 for fully developed 
turbulence in fluids). Energy in this terminology strictly refers to the energy in the Fourier spectrum of the data. The 
scaling index is often calculated from a linear regression of the E (k) plot over a chosen linear range of wave numbers 
(see Section 2.2 for examples applied to solar active regions, where Abramenko (2005a) and Hewett et al. (2008) use 
3—10 Mm, see Figure 20). More power at small k (hence large spatial scales) results in a larger scaling index, and so 
large j3 is suggestive of increased complexity in the system. Georgoulis (2012) studied a sample comprising hundreds 
of solar active regions and showed many of them follow non-Kolmogorov power-spectrum scaling, with j3 >5/3. 
Extended to a multi scale approach, this method can be used to eliminate any background non-SOC component Hewett 
et al. (2008). Fractals are defined in a similar manner as the self-similarity of an image across all scale sizes, or the 
scaling index of any length, I, to area. A, 

A^l^ . (44) 

The fractal dimension, a, and various other forms of fractal dimension (see McAteer (2013) for a complete list), is 
often calculated via a thresholding and contouring approach. The more complex the thresholded contour, the more 
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k (Mm"') 

Fig. 20.— The Fourier spectrum of a 2D slice of the active region magnetic field, plotted in log E (k) - log (k) show 
a clear linear range as a signature of SOC, and changes over shallow (27-Oct) to steep (30-Oct). From Hewett et al. 
(2008). 


space it fills, and therefore the larger the fractal dimension. McAteer et al. (2005) and Conlon et al. (2008) use such 
an approach to study the complexity of solar active regions. Georgoulis et al. (2002) adopt a similar approach to show 
the dust-like nature of small scale brightenings. Kestener et al. (2010) and Conlon et al. (2010) extended this to a 
multifractal approach that can be used to eliminate non-SOC backgrounds from images to show a clear relationship 
between the remaining multifractal spectrum of an active region and its potential to produce large solar flares. The 
power of these approaches, as evident in Figure 7 and Figure 20 is that they may provide a means of linking the clear 
time-varying nature of SOC avalanches in the emission from an active region (McAteer et al. 2007; McAteer and 
Bloomfield 2013) with a 2d spatial slice of the 3D SOC nature of spatial structures. However, it is important to note 
that although turbulence and fractality may be a signature of an SOC system there may be several other reasons for 
their occurrence. Therefore, these techniques should be accompanied by studies in time to confirm the existence of 
SOC (McAteer 2015). 


3.2. Feature Detection in the Temporal Domain 

An SOC system inevitably results in a series of catastrophies or avalanches, detectable in both observational data 
and simulations a release of energy. In an idealized dataset, each event would be well separated in space and time. 
A scientist simply needs to only identify each event, and can be secure in the knowledge that there is no overlap. 
However, such an idealized dataset is rare. Instead data often contains events that overlap significantly. In such a 
case of pulse-pile up, it may still be possible to separate out the signature of each individual event, and study these to 
determine if the waiting time distributions are the unique signature of SOC, or otherwise. 
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3.2.1. Power laws 

The Fourier power spectrum is a useful and simple tool to examine event occurrence in the temporal domain. 
Many systems exhibit power spectra such that the power spectral density P{v) is proportional to a negative power law 
of frequency V, 

p{v) oc v-P (45) 

where p >0. A nomenclature for noise spectra has emerged depending on the value of the index p, and is described 
in Table 2. Flicker, or shot noise, is common in electrical signals, and it was the analysis of this noise that produced a 
physically based model that is highly relevant for SOC models. Briefly, we envisage the electrical signal in an RCL- 

Table 2: Nomenclature of noise spectra (Aschwanden 2011) 


Index of power-law 

P 

Spectrum 

Nomenclature 

0 

white noise 

1 

pink noise, shot noise, flicker noise, 1 //noise 

2 

red noise, Brown(ian) noise 

3 

black noise 


circuit as consisting of the superposition of current spikes, parameterized as Dirac 5-functions having random arrival 
times tj, i.e., 

= (46) 

./ 

The general autocorrelation function given in Eq. (1) can be rewritten for a such a time series /(f) as 

C(f')= lim - / I{t)I{t + t')dt. (47) 

T j-T/i 

The Weiner-Khinchin theorem (Chatfield 1996) states that the power spectra density P(v) of a stationary random 
process is the Fourier transform of the corresponding autocorrelation function, 

/ -(-CXD 

dt'. (48) 

This enables the calculation of the power spectra of models of random processes /(f). Ziel (1950) and Aschwanden 
(2011) use the current model above to derive Schottky’s result (Schottky 1918) for the white noise spectral power 
distribution in electrical circuits. This general procedure in going from a model of the process to its power spectrum 
is used below to generate other power-law power spectra. 

Power-law power spectra have been observed in solar phenomena. McAteer et al. (2007) find power laws in solar 
flare X-ray data, and they then use this a means of studying the source of these X-rays in McAteer and Bloomfield 
(2013). Auchere et al. (2014) observe power laws in the integrated emission of small portions of active regions and 
the quiet Sun as observed in the 195A passband images from EIT over the frequency range 0.01 - 1 mHz. Ireland 
et al. (2015) observe power laws in power spectra of AIA 171A and AIA 193 A in active region, moss and quiet Sun 
areas in the frequency range 0.5 - 10 mHz. Gupta (2014) showed power-law power spectra in the intensity at six 
single points in AIA 171A coronal plumes extending over the frequency range 0.3 —4.0 mHz. Eurther out in the 
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solar atmosphere at 2.1 Rsun, Bemporad et al. (2008) show the presence of power-law power spectra in Ultraviolet 
Coronagraph Spectrometer observations of the intensity of Lyman-a in the frequency range 2.6 x 10^^ —1.3 x 10^^ 
Hz. Lower in the solar atmosphere, Reardon et al. (2008) show the presence of power-law Fourier power spectra, in 
the range 7-20 mHz, in the Doppler velocity of the chromospheric Ca II 854.2 nm line. 

Many models can generate power spectra that exhibit power laws. One simple model is the autoregressive process. 


Xt = aX,^i+N{0,a) 


(49) 


for f > 1, a > 0 and Gaussian noise A^(0, O’) with zero mean and standard deviation a. This simple process generates 
a power spectrum with index p = 2 in the limit of high frequencies (Chatfield 1996). Aschwanden (2011) gives 
the example of a shot noise spectrum of exponentially decaying pulses. Each pulse is modeled as an exponentially 
decaying function of time t 

f{t) = ^e^T, (50) 

for some timescale T and energy E. The corresponding Fourier power spectrum is (using Equations 47 and 48) 


l + {2nvT)^' 


(51) 


The total Fourier power spectrum of a distribution N{T) of these decaying pulses is 


P,o,ai{v)=Y^N{T)PT{v). (52) 

T 

Further, if the number of events of a given energy E is assumed to be 

N{E) oc (53) 

and the total energy in each event depends on its time scale T such that 

E oc T^+y (54) 

then it can be shown that the observed power spectrum can be approximated by 

Erom/(v)oc v-P-“£){i+r). (55) 

This derivation shows it is possible to generate a power law using swarms of statistically similar events. 

A power law may seem evident from a simple plot of the data, but the determination that a power law is actually 
present in the data is a subject that requires much attention (Clauset et al. 2009). There are essentially two parts to 
determining the properties of a power law in the data. Firstly, one must determine that a power law is an appropriate 
representation of the data. This should involve a combination of testing many different parameterizations/models of 
the data and making some determination as to which one best explains the data. Choosing a model also requires that 
the researcher think about the physical processes that may be occurring to generate the observations (Parnell 2002; 
Newman 2005; Vaughan 2010). This could be roughly classed as the model selection stage. The second stage is 
to actually determine the values of the parameters of the power law, properly taking in to account the details of the 
instrument, observational effects, and the statistics of the measurement. This is the parameter estimation stage. Clearly 
the two stages are intertwined to some extent. 

The identification that a power law is better than other reasonable models for the data is discussed by Clauset 
et al. (2009). A summarized procedure for deciding if a given data set follows a power law is given. This procedure is 
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applied to twenty four real-world datasets, drawn from a broad variety of disciplines, including physics, earth sciences, 
biology, ecology, paleontology, computer and information sciences, engineering, and the social sciences to test for the 
presence of power laws. The paper finds that in general, it is extremely difficult to tell the difference in the data 
between a log-normal behavior and that of a power law. 


Estimation of the parameters of a power law, along with an error estimation, is often crucially important as the 
index of the power law is often used as an indicator of the underlying physical process. Fitting a straight line to binned 
data is not recommended, as it introduces an arbitrary parameter the histogram binsize. White et al. (2008) show that 
no histogram binning yields values of the index of the power law consistently close to the true value. However, better 
methods exist. Let us assume a set of observations X, 1 < f < A^, drawn from a power-law probability density function 
of the form p{X) Newman (2005) derives that the likelihood function for this set of observations is given by 


7=1 + 


1 


N ^i=l 


InT,’ 


(56) 


There is an implicit assumption here that the observer has collected all the data perfectly, which is rarely the case. It 
can be appropriate to more carefully consider how the observation is made, and to include that in the estimation of the 
index of the power law. For example, Parnell and Jupp (2000) consider the observation process in the determination 
of the distribution of small heating events in the solar corona. It is assumed that the observed energy E^tis of a small 
heating event is related to its true energy E by 

Eobs = uE (57) 


where u is an under-reporting factor which satisfies 0 < u < 1. After some assumptions on the distribution of u, the 
under-reporting of the true energy of the event can be compensated for in the analysis; the final value of the index 
of the power law fully incorporates the modeling of the under-reporting. The likely physical nature of the energy 
deposition has also been considered for the same problem of the distribution of the energy in heating events in the 
corona. McIntosh and Charbonneau (2001) consider the geometry of the energy deposition event in the corona, which 
is shown to have a strong influence on the final value of the index of the power law. These two studies show that a 
careful consideration of the likely physical process, and the way it is observed, is required in order to fully realize the 
potential of the data. 


3.2.2. Pulse-pileup effects 

One of the tenets of slowly-driven SOC models is the separation of time scales, which means that the waiting 
time (i.e., the time interval between the starting times of two subsequent events) is larger than the event duration of 
the first event, so that there occurs only one event at a time, while no two events overlap with each other. While this 
requirement can easily be controlled in numerical cellular automaton simulations, it cannot be taken for granted when 
an automated pulse detection algorithm is applied to a time series of observations. In principle, numerical detection 
schemes can be designed to end one event before the next is detected, but this may truncate the duration of the earlier 
event or ignore a later event that starts during the decay phase of the earlier event. In practice, it is expected that 
the time separation criterion will be fulfilled during quiescent periods with low event rates, but it is possible that 
events start to overlap during more active periods, an effect known as pulse pile-up. This effect can be investigated by 
considering solar flare statistics during various phases of solar activity. 

Aschwanden and Freeland (2012) studied flare statistics from the GOES satellite sampled over a period of 37 
years (1975-2011), covering about three solar cycles. The soft X-ray flux from the Sun varies by about two orders of 
magnitude during each solar cycle, due to the variation of emerging magnetic fields and the resulting coronal plasma 
heating rate, which is all driven by the solar magnetic dynamo. This makes the Sun an ideal system to study SOC 
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Fig. 21.— Variation of the index of the power law, «/>(?), of the soft X-ray 1-8 A peak flux (top panel) and the flare 
rise time o:r(f), detected with GOES (middle panel), and the annual variation of the number of flares over 3 solar 
cycles (bottom panel). The flare rate predicts the variation in arit) of the flare time duration (smooth curve in middle 
panel) as a consequence of the violation of the separation of time scales. From Aschwanden and Freeland (2012). 


systems with variable drivers. While the power law of the soft X-ray peak rate is invariant, = 1.98 ±0.11 (Figure 21 
top), during different solar cycles, the time durations do have a variable slope from aj ~ 2.0 during solar minima to 
aj- « 2 — 5 during solar maxima. This is explained in terms of a flare pile-up effect. The variability of the flaring rate 
is shown in Figure 21 (bottom), from which the steepening of the index of the power law can be estimated by using 
the ratio of the mean inter-flare time interval to the mean flare duration (Figure 21 middle panel, solid curve), agreeing 
with the variability of the observed flare rate (Figure 21, middle panel, histogram). Apparently, the long flare durations 
are underestimated due to subsequent flares that start during the decay phase. This also affects the statistics of waiting 
times accordingly. In other words, the separation of time scales (i.e., the waiting times and flare durations) is violated 
during the busy periods of the solar cycle maximum. 

The influence of different pulse detection methods on the shape and index of the power law has also been studied 
in Buchlin et al. (2005), who compares a peak detection method, a threshold method, and a wavelet method. The 
peak method requires a relatively noise-free smoothed time profile, so that noise fluctuations do not contaminate the 
statistics with multiple peaks per time structure, leading to an excess of short waiting times. The threshold method 
requires that the time profiles return to a sub-threshold background level for each event, otherwise events in the 
decaying tail of a pulse time profile are ignored. The wavelet method has the ability to detect simultaneous pulses 
with different time scales, which would be impossible with the peak or threshold method. Interestingly, the three 
methods reveal quite different waiting-time distributions in each case. The threshold-based method seems to produce 
distributions that resemble power laws, while the peak-based and wavelet-based methods produce exponential-like 
distributions, at least in the regime of large waiting times. This result imposes some ambiguity in the interpretation of 
waiting-time distributions. The effect of event definition on the distribution of waiting times has also been numerically 
simulated with the continuously driven Olami-Feder-Christensen (OFC) model (Olami et al. 1992) by Hamon et al. 
( 2002 ). 





3.2.3. Waiting-time distributions 


In cases where pulse pile up can be neglected, or at least estimated and removed, it is possible to then study 
the waiting times between events as a possible signature of SOC. This leads naturally to the following key questions: 
do waiting-time distributions (WTDs) comprise an indisputable SOC-state feature? Can physical systems exhibiting 
different WTDs from the ones predicted in the original SOC concept be safely excluded from the long list of potential 
SOC systems? Since the development of the first avalanche models, it was suggested that the associated exponential- 
function WTDs should convey a necessary SOC signature. The context of solar flare dynamics provides a useful insight 
into this debate. Numerous researchers analyzed hard X-ray flare data in an attempt to construct the corresponding 
WTDs. Their results were initially conflicting. Biesecker (1994) used 1 yr of Gamma Ray Observatory (GRO) 
BATSE data to produce a WTD. The observed distribution was essentially exponential, covering the gaps due to lack 
of observational data through a simulation representing a Poisson process with a time-varying rate. Pearce et al. 
(1993), however, using 10 yr of Solar Maximum Mission hard X-ray burst spectrometer (HXRBS) data, found a WTD 
that was closer to a power law than to an exponential. This result suggested that the HXRBS events are interdependent. 
Crosby (1996) reported a distribution over a wide range of waiting times that could be fitted by a power law with an 
exponential rollover based on hard X-ray events observed in a single active region by the WATCH experiment onboard 
the GRANAT satellite. The index of the power law was close to that found by Pearce et al. (1993). 

Faced with these apparently conflicting results, Wheatland et al. (1998) re-examined the WTD of solar flare hard 
X-ray bursts. The WTD constructed from the ICE/ISEE 3 data showed an overabundance of short waiting times (10 s 
- 10 min) in comparison to a simulation of the time history of bursts as a Poisson process. This over-clustering with 
respect to a Poisson process indicates, according to Wheatland et al. (1998), the interdependence of some of the bursts 
that occurred in temporal proximity. Such a Poisson process would yield an exponential distribution for the waiting 
times of the solar flares and, according to Boffetta et al. (1999), such a distribution would only be expected if the events 
were completely uncorrelated. Moreover, Boffetta et al. (1999) suggested that SOC models are expected to display an 
exponential WTD P{Xl) =< Tl > ' < Tl >)> where < Tl >, is the average waiting time, which depends on 

the parameters of the model. This behavior is related to the fact that the avalanche duration is much smaller than the 
loading time {i.e., the time between two successive injections of magnetic field in random positions) and charging place 
(i.e., the random position in which the injection of the magnetic field takes place) is independent from the avalanche 
position. Then one expects no correlation between successive bursts and thus a trivial, exponential statistics for the 
waiting times. However, various caveats on this assessment were thereafter voiced: first, Buchlin (2005) suggested 
that thresholding the event time series may result in WTD resembling power laws in an SOC system. Based on a 
non-stationary Poisson model as introduced by Wheatland (2000) and further discussed in Wheatland and Litvinenko 
(2002), Aschwanden and McTiernan (2010) reviewed numerous studies and data sets to conclude that WTD for solar 
flares can generally be approximated by a non-stationary Poisson distribution of the form P(At) Ao(l -f AqA?)^^, 
where Ao = 1/Afo is the flare rate corresponding to a waiting time Afo, below which there is a high flare rate, or 
clustering, of small released energies. Above this time, the flare rate decreases with flare magnitude (released energy) 
giving rise to a WTD that resembles a power law. Evidence that this WTD can in fact correspond to an SOC system 
also stems from the analytical predictions of the avalanche model of solar flares Charbonneau et al. (2001) and the 
fractal-diffusive SOC model described byAschwanden (2014) and discussed extensively by Aschwanden et al. (2014) 
(this volume). 

Nonetheless, Boffetta et al. (1999) calculated the waiting times for flares recorded in hard X-rays during the 
period 1976-1996. Two different datasets were created: dataset A, by calculating only the differences between the 
time of occurrence of flares within the same active region and dataset B, by calculating the time differences between 
two successive maxima of flare intensity regardless of the position of the flare on the Sun’s surface. The results 
presented in Figure 22 distinctively show a power-law distribution of WTDs for both datasets A and B. In the inset 
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of this figure Boffetta et al. (1999) show the WTD distribution for dataset B (solid line) derived from observations, 
compared with the corresponding distribution obtained through a cellular automaton model (dashed line) used by 
Boffetta et al. (1999) as reference of the exponential behavior of SOC simulations. 

These results have heed used to argue against the relevance of SOC in solar-flare dynamics. It has been also 
proposed that SOC should be discarded in plasma turbulent transport dynamics in magnetic confinement devices 
after carrying out the same analysis on edge electrostatic fluctuations from the reversed-fleld experiment (RFX) pinch 
(Spada et al. 2001). Yet, as suggested by Sanchez et al. (2002), such tests must be considered with extreme care. In 
their work, Sanchez et al. (2002) stressed that the waiting time definition is of crucial importance with regards to the 
resulting WTD of a physical or simulated system. Until then, some authors used the time interval between triggers, 
others the time interval between two consecutive maxima in burst intensity, and Anally others considered the time lapse 
between the end of a burst and the beginning of the next one. Sanchez et al. (2002) showed that only the quiet time 
would yield an exponential WTD for non-correlated triggers in an SOC system. Sanchez et al. (2002) carried out their 
simulations on a ID running sandpile, consisting of L cells, and with a closed and an open boundary, respectively. 
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Fig. 22.— Probability distribution function of the waiting time P{'Cl) between two X-ray flares for two datasets (A, 
dashed line and B, solid line). The straight lines are the respective least-squares fit of a power law. The inset shows 
the distribution for dataset B (solid line) and the distribution obtained through a reference SOC model (dashed line) 
that exhibits an exponential distribution. The variables shown in the inset have been normalized to the respective 
root-mean-square values. From Boffetta ef aZ. (1999). 
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located at the first and last cells. At each iteration, Uq grains of sand are dropped at each cell with probability Pq. 
Whenever the local sand slope, Zj = hj — exceeds some prescribed critical value Z^, Nf grains of sand are moved 
to the next cell. The sandpile reaches the critical state after the incoming sand flux is balanced by the flux leaving the 
system through the open boundary. With these results, Sanchez et al. (2002) claimed that the lack of an exponential 
WTD should not be used to discard SOC dynamics when all other signatures (i.e., regions in fluctuation power 
spectra, or Hurst exponents H > 0.5 from the rescaled range (R/S) analysis) suggest the existence of SOC. They 
propose that an exponential WTD is not a necessary condition for SOC state in the following cases: 

a) When the avalanche durations are longer than the quiet times; then power laws can appear because waiting times 
become contaminated by the event-duration scaling go the power law. 

b) When the avalanche durations are much shorter than the quiet times; then power laws can still appear if the 
measurements’ maximum resolution lies within the self-similar range, since all detected avalanches then become 
strongly correlated. This argument was supported by the earlier study of Christensen and Olami (1992), in the 
context of a spring-block model for earthquakes. The model showed that waiting times could follow power-law 
distributions in case events larger than a certain size are only considered. 

c) When experimental resolution is sufficiently high to detect events of all possible sizes; the lack of exponential 
waiting times in this case might simply imply that the system is driven in a correlated way. The physical origin 
of the correlated driver in this case is system-dependent and should be determined on a case-by-case basis. 

It is therefore possible that a system governed by SOC dynamics can lack exponential WTD statistics, not only when 
the experimental resolution lies within the self-similar scale range, but also when the system is slowly driven in a 
correlated way. Appreciating the long-standing debate at this point, we recommend caution in the interpretation of a 
given WTD and suggest that waiting-time statistics should not be used as a necessary test of SOC behavior in physical 
systems. 


3.3. Feature detection in the Spatial-Temporal Domain 

The previous two sections have focused on identifying features either in space or in time. This is appropriate as 
scientists are often relegated to studying such datasets. A time series is often all that is obtained from stellar observa¬ 
tions. Although this can reveal time-separable pulses that can be used for testing the statistics of SOC phenomena, all 
spatial information is concealed in a dot-like point source. More informative from imaging observations can exhibit 
the detailed fractal spatial structure of SOC phenomenon, but temporal information is commonly lacking or ignored. 
Combining the two domains of space and time into spatio-temporal event detection methods clearly present a powerful 
means to analyze SOC phenomena. However, these methods are quite complicated and hence need a sophisticated 
initial setup in order to work correctly. 


3.3.1. Spreading and Avalanche Exponents 

The relationships between the spreading and avalanche exponents (Munoz et al. 1999) and spatio-temporal struc¬ 
tures provides a useful method to study if the system is in an SOC state. The concepts of spreading and avalanche 
exponents were put to use in the case of numerical models for magnetospheric (Morales and Charbonneau 2008a) and 
solar flare (Morales and Charbonneau 2008b) phenomena as well as with observations of auroral emissions (Uritsky 
et al., 2000) and multi-wavelength data for solar flares (Aschwanden 2012). When an SOC system arrives in the 
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vicinity of criticality the spreading of an active site can be described by a number of scaling laws that characterize its 
dynamical properties. Generally the measured quantity is a survival probability P{t) that an instability is still active 
after t iterations and the number of active sites at a given time, n{t) (Bonachela and Munoz 2007). Both quantities are 
expected to satisfy a power law with f, 

(58) 

where 77 and d are the so-called spreading exponents (Munoz et al. 1999). This implies that the total number of active 
sites having a lifetime T scales as n* ^ jn+5 ^ therefore its time integral should be characterized by the exponent 
K= 1+77 + 5 . Provided that these scaling relations hold, then the total number of avalanching sites - the size of the 
event - 5, scales with its lifetime T as: 

S{T) - . (59) 

Another spreading exponent that characterizes the probability distribution of avalanche sizes is therefore found as 

As avalanches of size S can have different durations T, the probability of an avalanche reaching a size s before 
dying is 

r^max ^ 

P{s)= P{s\t){l-t-^)dt , (60) 

^ ^min 

where and t^ax are the upper and lower duration bounds of size-s avalanches, and P(s|7) is the conditional proba¬ 
bility of an avalanche having reached size s at time t since onset. P(s|f) is bell-shaped and peaks at 7 ~ l/s*+')+^ so it 



Fig. 23.— Correlation plot of avalanche sizes (S) vs lifetimes {T) for a simulation on a square lattice of size A = 128 
and angular threshold 0^ = 2.25 rad. The gray line is a least-squares fit, computed using only avalanches with lifetime 
r > 40 iterations. The value of the spreading exponent in this case is if = 1.82 + 0.3. As 5 = ~ T'^,L ^ 

which is this case is L ~ T’O.ei^ close to classical diffusion (L ~ From Morales and Charbonneau (2008a). 
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can be shown (Munoz et al. 1999) that P{s) scales as 

P{s) , 

with the same scaling as expected for P{S). 

These redundant relations provided in Equation 60 and 61 can be considered as another way of verifying if an 
avalanching system is in an SOC state. These relations were confirmed and presented for the case of an anisotropic 
SOC model for solar flares that used magnetic field lines as a basic dynamical element, and the angle between field 
lines as the threshold value (Morales and Charbonneau 2008b). The typical correlations found between the avalanche 
sizes and lifetimes are displayed in Figure 23. The same analysis has also proved useful for the case of an SOC model 
for the magnetosphere (Liu et al. 2010). In the last decades it has been claimed that the solar corona and the Earth’s 
magnetosphere might be in SOC. Several models have been produced in order to prove this assertion and the formalism 
of spreading exponents indeed provides an excellent venue to test observational data and models. 




l + ri + 25 
1 + 77 + 5 


(61) 


3.3.2. Spatio-Temporal Structures 

Spatio-temporal structures are well defined in classical SOC models, such as a numerical cellular automaton 
simulation like the BTW model (Bak et al. 1987). Once an SOC avalanche starts at time fi, the evolution of the 
avalanche size is updated as described in Section 2.3 above. In this section we describe the spatial-temporal evolution 
that determines the resulting size off the avalanche.The initial size of the avalanche at time fi has then the size s, = 1, 
which represents the unstable node in the lattice grid. In the next time step, zero to four next neighbors can become 
unstable (in a 2D lattice grid), after the application of the SOC re-distribution rule, and thus the avalanche has a size of 
S 2 = 1, • ■ •,4 nodes, or dies out (^2 = 0). If the avalanche is further unstable, the size can grow to s, = 1,2,..., 8(/ >2) 
next neighbors, and so forth. The cumulative avalanche size after time step is the time-integrated instantaneous size 
of the avalanche, i.e., 

S= / s{t)dt = ysi. (62) 

Jo 

If the same spatial pixel is active multiple times during an avalanche event, it is counted multiple times correspond¬ 
ingly. Consequently, the so-defined avalanche size S is not a geometric volume, but rather a volume in hyper space 
(with d geometric dimensions plus one time dimension). Note that the time step and the spatial pixel (or voxel) size are 
dimensionless in numerical lattice simulations and are set to unity for convenience. Non-imaging astrophysical obser¬ 
vations typically record the spatio-temporal information of an SOC phenomenon by a flux or intensity 7) = F{t = 7,) 
at time 7, with a cadence or time interval dt. The summed flux adds up to a time-integrated fluence or energy E as 
discussed in Section 3.1.2. The flux Fi corresponds to the emission from all active or unstable pixels in a cellular au¬ 
tomaton avalanche, and thus represents the instantaneous energy dissipation rate dEt/dt at time 7,. The total dissipated 
energy per avalanche, E, corresponds then to the time-integrated size S, as E — AE S S, with a constant energy 
dissipation quantum AE per pixel or voxel. 

With the luxury of high-resolution imaging when observing a candidate SOC phenomena in astrophysical data, 
it is possible to additionally measure the (possibly fractal) area A, at each time 7,. This renders a snapshot of the 
instantaneous contours of an SOC avalanche, defined by the sum of pixels with a flux in excess of some noise threshold, 
Fi > Ffh (similar to Sections 3.1.2 and 3.1.3). The area information A (7) is not sufficient to reconstruct the volume 
y(7) at a given time 7, as there is no direct information on the column depth along the line-of-sight. However, the 
information of the avalanche location is crucial to separate multiple avalanches occurring at the same time, or over¬ 
lapping in time, at different spatial locations. The concept of the spatio-temporal tracking of two time-overlapping 
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Fig. 24.— left: An example of a POLAR UVI image. Right: A schematic drawing illustrating the method of identi¬ 
fying spatio-temporal auroral events from POLAR UVI images. The elliptical spots in the image planes indicate the 
time evolution of two time-overlapping auroral events with the photon flux exceeding some noise threshold. From 
Uritsky et al. (2002). 


avalanches is shown in Fig. 24 for the case of auroral sizes recorded with POLAR UVI (Uritsky et al. 2002). The 
area distribution of auroral sizes was found to have a power-law distribution with a slope of Ua = 1.73 ±0.03 for the 
auroral observations during Jan 1997. In contrast, earlier measurements by Lui et al. (2000) of the same data yielded 
a much flatter distribution with a slope of = 1.21 ±0.08, because multiple time-overlapping auroral events were 
not spatially separated, and thus led to an over-estimation of large areas. The flatter index is also not consistent with 
predictions of a theoretical SOC model (see Section 3.3.1 in Aschwanden et al. (2014), this volume). Therefore, the 
proper spatial separation of time-overlapping events in spatio-temporal detection methods is very important to obtain 
the correct SOC statistics. 

Spatio-temporal detection of nanoflares in the solar corona present a good example of the power of this technique. 
Nanoflares often occur near-simultaneously in different spatial locations, and thus require a sophisticated automated 
feature detection algorithm. While an absolute flux threshold, i.e., 7) > 7)/,, was used in the foregoing description of 
detecting auroral events, solar nanoflares cannot be detected by an absolute flux threshold, because they are associated 
with much weaker and fainter local brightness enhancements than the variation of the flux in the surrounding or co- 
spatial active regions, or quiet Sun. Active regions might have a brightness of 7" « 10^ — 10® DN/s in typical EUV 
images, while nanoflares exhibit only tiny brightness variations in the order of 7" « 1 — 10^ DN/s (Aschwanden et al. 
2000a,b). Nanoflares therefore have to be detected by their temporal variability, rather than by their absolute flux: 
consequently a time variability threshold between two consecutive images should be applied. 


F{x,y\ti+i)-F{x,y\ti) > M^,hresh = 3o/ , (63) 

rather than an absolute flux threshold. A possible threshold, e.g., Ff^resh = 3(1/, can be specified by the photon Poisson 
noise in a time bin, with additional correction for spatial rebinning (to macropixels), exposure time, and other instru¬ 
mental effects. An example of a solar EUV image is shown in Eigure 25, where the detected nanoflares are marked 
with ellipses.The location of detected nanoflares are not necessarily coincident with the locations of highest bright- 
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Fig. 25.— Spatial maps of 20 EUV nanoflare events are shown, observed with TRACE in 195 A on 1999 Eebruary 
17, 02; 16-02:59 UT. The greyscale images (first and third column represent difference images taken at the peak and 
minimum time of each nanoflare and averaged over five cadences. The contours of these difference images of detected 
nanoflares (second and fourth column) have a flux increment of 4 DN. Erom Aschwanden et al. (2000b). 


ness, but their flux variability exceeds a threshold F > Fthresh in a difference image. Examples of variability maps 
are shown in Figure 26 which show the contours of EUV brightness, the pixels with significant variability (crosses). 
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195 A, event # 0 



195 A, event # 3 



195 A, event # 6 
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195 A, event # 9 



195 A, event# 1 


195 A, event #2 



Fig. 26.— Spatial clustering of the pattern recognition code is illustrated for the 12 largest events on 1999 February 
17, 02:15-03:00 UT. The contours outline local EUV intensity maps around the detected structures. The crosses mark 
the positions of macropixels with significant variability (A^ < 3(7). The spatiotemporal pattern algorithm starts at the 
pixel with the largest variability, which is located at the center of each field of view, and clusters nearest neighbors 
if they fulfill the time coincidence criterion. These macropixels that fulfill the time coincidence criterion define an 
event, marked with diamonds, and encircled with an ellipse. Each macropixel that is part of an event, is excluded in 
subsequent events. Erom Aschwanden et al. (2000a). 


and pixels with significant variability that is cospatial in two subsequent images (diamonds). The automated detection 
criterion needs to include both spatial coherence and temporal contiguity. Those pixels that fulfill both criteria are 
marked with an elliptical area A that characterizes the Euclidean flare area, while the diamonds in Eigure 26 demarcate 
the instantaneous fractal flare area. 
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The numerical event detection code used for the examples shown in Figure 25 and Figure 26 was especially 
designed to detect solar microflares and nanoflares, which represent the faintest counterparts of solar flares, and thus 
are important to extend the dynamic range of frequency distributions of flare energies over nine orders of magnitude. 
Similar codes were also developed by Krucker and Benz (1998) and Parnell and Jupp (2000), which triggered con¬ 
troversial results on the index of the power law in the nanoflare regime. A number of assumptions were considered 
that contribute to the initially discrepant results of these indices, such as event definition, selection, and discrimina¬ 
tion, sample completeness, observing cadence and exposure times, pattern recognition algorithms, threshold criteria, 
instrumental noise, wavelength coverage, fractal geometry, but also physical modeling issues of energy, temperature, 
electron density, line-of-sight integration, and fractal volume (e.g., Aschwanden and Parnell 2002; Benz and Krucker 
2002). 
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4. Summary And Conclusions 

In this review we have shown that the numerical detection of SOC is a research field onto itself. Although it 
remains difficult to state definitely that a system exists in a state of SOC based on feature detection alone, much 
progress has been made across all science fields that set out to attempt this feat. The basic studies of autocorrelations 
provides a powerful tool to determine if a system is in SOC. It can be used to determine if the particles in the system 
are spatially and temporally correlated in the appropriate manner, and is readily applicable to both simulations and 
experimental data. The structure function provides a complementary method using field increments, and provides an 
analytical connection to studies of SOC geometry. Future progress will surely consist of combining such methods with 
the more application-oriented methods such as marginal stability and statistical stationarity to high spatial resolution 
data. Even when such data is not available, block scaling provide a powerful technique to extract potential signatures 
of SOC. 

The problems associated with working with less-than-optimum data are discussed in detail in Section 3. The 
scientist is reduced to applying some thresholds, and usually does not have all measurements in full 4 dimensions (3 
space and 1 time). However, even with static 2D spatial slices, progress in this field has been made by adopting and 
adapting techniques of detecting power laws and fractals. Such features are undoubtedly ubiquitous in nature, and 
may well be a good signature of SOC systems. However we urge caution in adopting either of these as being a unique 
signature of SOC without further independent studies. In particular, the detection of power laws has undergone its 
own revolution in the past few years and powerful statistical tools are now freely and widely available for all scientists 
to use. Combined with a full understanding of instrumental effects of sub-sampling of the system, this opens up future 
studies in waiting time distributions as a signature of SOC, especially in those areas of study with long, homogenous, 
uninterrupted datasets. In terms of identifying features, it seems clear that the confidence in assigning the label of SOC 
to a system is much greater when we include as many datasets as possible, and as many dimensions as possible. In 
particular, if data can be used to move from units of DN (or counts per second) to units of energy (or energy per second) 
we will undoubtedly obtain a better measure of the energy release processes. It is these energy release processes that 
we then attempt to recognize. Probably the greatest untapped potential for the next 25 years lies in spatio-temporal 
studies. The concepts of spreading and avalanche exponents can be adopted for all future datasets. As hi-fidelity, 
multi-spectral data becomes more commonly available across all areas of science, perhaps the biggest obstacle to 
success is the risk of a lack of the interdisciplinary research avenues (such as the ISSI workshops), necessary to help 
us exploit each others’ data. Numerical methods will play a key role in the advancement of clearly and unambiguously 
detect SOC in data. Scientists must continue to explore and understand these methods as applied to each others’ data, 
as numerical methods will surely continue to provide a key link between simulations and experiments across all fields 
in scientific research. Advances in any field of research must spread across all of science. We must continue to seek 
to explore this interdisciplinary boundary over the next 25 years. 
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